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SHAPE FROM SHADING: A METHOD FOR OBTAINING 
THE SHAPE OF A SMOOTH OPAQUE OBJECT FROM ONE VIEW* 


Abstract 


A method will be described for finding the shape of a smooth 
opaque object from a monocular image, given a knowledge of the 
surface photometry, the position of the light-source and certain 
auxiliary information to resolve ambiguities. This method is 
complementary to the use of stereoscopy which relies on matching 
up sharp detail and will fail on smooth objects. Until now the 
image processing of single views has been restricted to objects 
which can meaningfully be considered two-dimensional or bounded 
by plane surfaces. 

It is possible to derive a first-order non-linear partial differ¬ 
ential equation in two unknowns relating the intensity at the image 
points to the shape of the object. This equation can be solved by 
means of an equivalent set of five ordinary differential equations. 

A curve traced out by solving this set of equations for one set of 
starting values is called a characteristic strip. Starting one of 
these strips from each point on some initial curve will produce the 
whole solution surface. The initial curves can usually be construct 
ed around so-called singular points. 

A number of applications of this method will be discussed including 
one to lunar topography and one to the scanning electron microscope. 
In both of these cases great simplifications occur in the equations. 
A note on polyhedra follows and a quantitative theory of facial make 
up is touched upon. 

An implementation of some of these ideas on the PDP-6 computer with 
its attached image-dissector camera at the Artificial Intelligence 
Laboratory will be described, and also a nose-recognition program. 


*This report reproduces a thesis of the same title submitted to 
the Department of Electrical Engineering, Massachusetts Institute 
of Technology, in partial fulfillment of the requirements for the 
degree of Doctor of Philosophy, June 1970. 
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J . INTRODUCT ION : 

J.J SHADING AS A MONOCULAR DEPTH CUE: 

Consider a smooth object known to have a uniform surface. An 
image of such an object will exhibit shading (gradations of 
reflected light intensity) which can be used to determine its 
shape/ given only a picture from a single viewpoint. This is 
not obvious since at each point in the image we know only the 
reflectivity at the corresponding object point. For some 
points (called singular points here) the reflectivity does 
uniquely determine the local normal, but for almost all 
points it does not. The shape of the surface cannot be found 
by local operations alone. 



Figure 2 : Definition of the incident (i), emittance (e) 
and phase angle (g). 
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For many surfaces the fraction of the incident tight which is 
scattered in a given direction is a smooth function of the 
angles involved. It is convenient to think of the situation 
as depending on three angles: the incident angte (between 
local normal and incident ray)/ the emittance (or emergent) 
angle (between local normal and emitted ray) and the phase 
angle (between incident and emitted rays). 

It can be shown that the shape can be obtained from the 
shading if we know the reflectivity function and the position 
of the light-source(s ) . The reflectivity and the gradient of 
the surface can be related by a non-linear first-order 
partial differential equation in two unknowns. The recipe 
for solving this equation is to set up an equivalent set of 
five ordinary differencial equations (three for the 
coordinates and two for the components of the gradient) and 
then to integrate these numerically along certain curved 
paths on the object called c ha racte r i s t ics [5]. For while we 
cannot determine the gradient locally, we can/ roughly 
speaking/ determine its component in one special direction. 
Then taking a small step in this direction/ we can repeat the 
process - the curve traced out on the object in this mannet 
is called a character ist ic. Its projection on the image 
plane will be referred to as the base characteristic. The 
shape of the visible surface of the object is thus given as a 
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sequence of coordinates on some such curves along its 
surface. 

An initial known curve on the object Is needed to start the 
solution. Such a curve can usually be constructed near the 
singular points mentioned earlier using the known 'local 
normal. The only additional Information needed is the 

distance to the singular point and whether the surface is 
convex or concave w.r.t. the observer at this point - such 
ambiguities arise in several other instances in the process 
of solution as will be seen. 



Figure 3 : Image of a sphere and a stereo-pair of the 

characteristic curves obtained from the shading. 

To solve the equat ions, the reflectivity as a function of the 

three angles must be known, as well as the geometry relating 

light-source, object and observer. Multiple or extended 
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light-sources increase the complexity of the solution 
algorithm presented. But all of this initially needed 
information can be deduced from the image if a calibration 
object of known shape is present in the same image. 
Furthermore/ incorrect assumptions about the reflectivity 
function and the position of the light-source(s) can lead to 
inconsistencies in the solution and it may be possible to 
utilize this information in the absence of a calibration 
obj ect. 

In practice it is found that if the object is at all complex, 
its image will be segmented by edges. Some of these are 
purely visual/ due to the occlusion of one surface by 
another/ others are angular edges (also called joints here) 
on a single object. Another kind of edge is the ambiguity 
edge. This is an edge which the characteristics cannot 

cross, indicating an ambiguity which cannot be resolved 
locally. One can solve inside each region bounded by these 
various edges, but some global or external knowledge is 
needed to match up the regions. In the case of an angular 
edge on the object one can integrate up to the edge and then 
use the known location of the edge as an initial curve for 
another region (provided one resolves the ambiguity present 
here, as on all initial curves). 
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A very similar situation also obtains when one bridges a 
shadow. Since one edge of the shadow and the position of 
the light source is known, we can trace along the rays 
grazing the edge until the corresponding image points fall on 
an illuminated region. Since we know the path of each ray, 
we can calculate the coordinates of the point where It 
impinges on the object. The edge of the shadow (which need 
not be on the same object) can now serve as an initial curve 
from which to continue the solution. 

A number of interesting applications of this method can be 
mentioned. The first of these concerns the scanning electron 
microscope (SEM) which produces images v/hich are particularly 
easy to interpret, since the intensity recorded is a function 
of the slope of the object at that point and is thus a form 
of shading (as opposed to optical and transmission electron 
microscopes which produce intensities which depend on 
thickness and optical or electron density). The geometry of 
the scanning electron microscope a l l ows several 
simplifications in the algorithm for determining shape from 
shading (e.g. there are no shadows). Because of the random 
access capability of the beam of this microscope it should be 
easy and useful to combine it with a small computer to obtain 
three-dimensional information directly. 
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Another important application lies in the determination of 
lunar topography. Here the special reflectivity function of 
the material in the maria of the moon allows a very great 
simplification of the equations used in the shape-from- 
shading algorithm. The equations in fact reduce to one 
integral which has to be evaluated along each of a family of 
predetermined straight lines in the image, making for high 
accuracy. This problem was first tackled for areas near the 
terminator (the dividing line between the illuminated and the 
un i l l uminat ed part of the moon's disk) by J . van Diggelen at 
the Astronomical Institute of the Netherlands in 7957 [2] and 

solved by T. Rlndfleisch at the Jet Propulsion Laboratory in 
1966 [4] and the method applied to several pictures returned 
by the Ranger spacecraft. This gave the first indication 
that the general solution discussed here might be possible. 

It should be pointed out that this method is complementary to 
the use of stereopsis, since the latter will match up sharp 
detail and edges while the shading information will determine 
the shape of the smooth portions of the surface. 

So far we have assumed that the surface is uniform in its 
photometric properties. Any non-uniformity will cause this 
algorithm to determine an incorrect shape. This is one of 
the uses of facial make-up; by darkening certain slopes they 
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can be made to appear steeper for example. In some cases 
surface-markings can be detected if they lead to 
discontinuities of the calculated shape. 

Judging by our wide use of monocular pictures (photographs or 
even paintings and woodcuts) of people and other smooth 
objects/ humans are good at Interpreting shading information. 
The short-comings of our method which are related to the 
shading information available can be expected to be found in 
human visual perception too. It will of course be difficult 
to decide whether the visual system actually determines the 
shape quantitatively or whether it uses the shading 
information in a very qualitative way only. A quantitative 
determination would involve operations more complicated than 
those used in edge-finding for example. Since the 
information is not local/ the surface-shape calculations 
cannot be carried out entirely in parallel. 


1.2 HISTORY OF THE PROBLEM: 

After formulating the image illumination equation as the 
basis of a method of finding shape from shading/ a literature 
search was performed to see if a solution had previously been 
obtained. The literature on perception has only a few 



Pose 14 


conjectures on the possibility of determinins shape from the 
monocular depth-cue of shading. Photogrammetry does not pay 
much attention to the reflectivity function / but only various 
integrals of it, measured by such devices as the integrating 
photometer. With few exceptions machine perception so far 
has been restricted to objects which can usefully be 
considered two-dimensional and objects bounded by planes 
(polyhedra). 

The one relevant research was found in the paper on lunar 
topography by T. Rindfleisch [4] which gives complete details 
of a solution obtained in the form of an integral in the 
special case of the reflectivity function of the moon. This 
raised the hope that a general solution existed. The 
(x',y',r) coordinate system used in [4] leads to Intractable 
equations - but we found a solution using a different 
coordinate system, (x*,y*,z). As a check the solution for 
lunar topography was rederived from this set of equations 
(Rindfleisch found his solution in quite a different manner - 
searching for predetermined curves in the image along which 
the surface can be found as some integral involving the 
measured image illumination). A first program (old SHADE) 
v/as then written which solved along one characteristic at a 
time using various p red ict or-correct or-modi f ie r methods 17 3 • 
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Another program (RFFLEC) was used to measure the reflectivity 
function from a calibration sphere. Various short-conings of 
our image-dissector sensing device were affecting the 
accuracy of these measurements. Since very little was known 
about the characteristics of this device on other than 
theoretical grounds [9], a program (TFXTUR) was developed to 
measure various properties such as resolution, signal to 
noise ratio, drift, settling time, scatter and pinholes in 
the photocathode. An attempt was then made to provide 
software to compensate for some defects such as distortion 
and non-uniform sensitivity, using measurements from test 
patterns (DISTOR). 

These techniques allowed an estimation of what accuracy can 
be achieved under optimal conditions. The program had 
numerous problems when dealing with objects other than simple 
convex ones (mostly because it solved each characterIstic 
separately) and as must be apparent, was sensitive to 
imperfections in the sensing device (partly because of the 
way it obtained intensity gradients). 

After the defects in the first program had been found, and a 
decision made to rewrite it, a great simplification of the 
main equations was found using a different coordinate system 
(x,y,z) and a slight extension of standard vector notation 
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(the voluminous equations for the inconvenient coordinate 
system Cx',y',z) are not reproduced here). An unfortunate 
but unimportant side-effect is an increase in the complexity 
of the derivation of the lunar topography integral. The new 
equations and numerous changes in the method of solution were 
incorporated in a new program (new SHADE) which was less 
sensitive to the various shortcomings of our image-dissec tor. 
This program can handle objects somewhat more complicated 
than its predecessor and solves all characteristics at the 
same time. 

In parallel with the programming work, theoretical efforts 
were made to define and get around some of the difficulties 
of the method of shape from shading. Of particular interest 
were applications where the equations simplify greatly. 
Unfortunately the massive simplification found in the case of 
lunar topography is unique. Of most interest are the cases 
where we have some advance knowledge of the characteristics 
(for lunar topography they are completely independent of the 
image - for the scanning electron microscope they are paths 
of steepest descent). 



1.3 PREVIEW OF CHAPTERS - GUIDE TO THE HURRIED READER: 


References to articles and hooks listed at the end will be by 
numbers enclosed in brackets. Numbers contained in 
parentheses refer to sections and subsections in this work. 
In an attempt to be complete/ a few subsections were included 
which will have only limited appeal to some readers; hence 
this guide. 

Chapter 1 provides an introduction to the depth-cue of 
shading/ its use in determining shape and its history. 
Chapter 2 develops the necessary equations in detail/ 
starting with the definition of the reflectivity function. 
Subsections 2.1.2 to 2.1.4 and 2.2 can we l l be skipped by the 
hurried reader. In section 2.3 the partial differential 
equation is obtained, the vector di f fe rent iat i on notation 
introduced and an equivalent set of five ordinary 
differential equations derived. Section 2.3 is perhaps the 
most important section. Sections 2.12 to 2.16 deal with some 
miscellaneous implications and may be omitted without loss of 
contin uit y. 

Chapter 3 describes in detail some practical situations where 
the special conditions encountered make use of the method of 
determining shape-f rom-shad ing particularly attractive. 
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Section 3.1 deals with the scanning electron microscope. The 
reader should be warned about the tedious derivation of the 
simple integral for the case of lunar topography in 3.2 . 
Omitting subsection 3.2.3 will avoid the bulk of the 
algebraic detail, and most of the conclusions will be found 
anyway in subsection 3.2.4 . 

Chapter 4 describes the experiments carried out with the two 
programs (using the results developed in chapter 2) to obtain 
shapes from images projected on an image dissector camera 
attached to the PDP-6 computer in the Artificial Intelligence 
Laboratory. Section 4.1 deals with the less successful first 
program, and contains details on auxiliary routines. Section 
4.2 deals with the second program wh i ch solves the 
characteristics in parallel and also uses the important 
sharpening process. Sections 4.1 and 4.2 are next in 
importance to section 2.3 . 

Section 4.3 describes an application to a recognition task - 
that of nose-recognition. Section 4.4 contains an overall 
summary and conclusions about the capabilities of the method 
of shape-from-shading, with subsection 4.4,1 giving 
suggestions for future investigations. This is followed by a 


list of references. 
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2. THEORETICAL RESULTS: 


2 .1 THE REFLECTIVITY FUNCTION: 

2 .1.1 DEFINITION OF THE REFLECTIVITY FUNCTION: 



Figure 4 : Illustration of the variables used in the 
definition of the reflectivity function. 

Consider a surface element of size dS inclined i w.r.t. the 
incident ray and e w.r.t. the emitted ray (The angles are 
measured w.r.t. the normal). Let the incident light 
intensity be I, per unit area perpendicular to the incident 
ray. The amount of light falling on the surface element is 
t hen Ij cos(i) dS. 

Let the emitted ray have intensity I * per unit solid angle 
per unit area perpendicular to the emitted ray. So the 
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amount of light intercepted by an area subtending a solid 
angle dv. at t he surface element will be I z cos(e) dS dw. The 
reflectivity funcion 4>(i,e,g) is then defined to be I x /I,. 

If we want to be more precise about what units the intensity 
is measured in, we have to take into account the spectral 
distribution of the light emitted by the source/ as well as 
the spectral sensitivity of the sensor (with this proviso we 
can speak of watts per unit area and watts per unit solid 
angle per unit area etc.). We need not be too concerned with 
this if we either use white paint, or measure the 

reflectiivity function with the same equipment later used in 
the shape-from-shading algorithm. It should be noted that 
for most surfaces the reflectivity function is not 

independent of the color of the light used. Typically the 
specular component of the reflected light, being reflected 
before it has penetrated far into the surface, will be 
unchanged, while the matt component will be colored by 
pigments in the surface coating. 

Several other definitions of the reflectivity function are in 
use which are multiples of the one defined here by TV, 2, 
cos(e) and/or cos(i). The specific formulation chosen here 
makes the equation relating the incident light intensity to 
the image illumination very simple. 
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2 . 1.2 FUNCTIONS DERIVED FROM THE REFLECTIVITY FUNCTION: 

The next few subsections (2 .1 .2 .1 , 2 . 1 . 2.2 and 2 . 1 . 2 . 3 ) are 

Included to relate the reflectivity functions to those more 
commonly mentioned in the literature. Some readers may want 
to skip these subsections. 

2 . 1 . 2.1 THE INTEGRATING PHOTOMETER: 

A flat sample of the surface under investigation is mounted 
in the center of a hollow sphere coated on the inside with a 
highly reflective matt substance. Through one small hole a 
light ray enters and impinges on the sample with incident 
angle i. A photosensitive device is introduced through 
another small hole and measures an intensity proportional to 
the light scattered by the sample into all directions. 

The total intensity measured is: 




lonio-photometer (used for measuring 
tivity functions). 



ec 


SENSOR 
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tv/?. 

I^cosCe) dS sin(e) de dA 
t r/o. 

c|>(i,e,g) J f cos(e) dS sin(e) de dA 

Where cos(g) - cos(i) cos(e) + sin(i) sinCe) cos (A) 

The total incident intensity Is I cos ( i) dS. The fraction of 
light reflected is then: 

y'-siTT >-xTT/2. 

b (i) « [ J I ( i / e, g )* (I / 2) * s i n ( 2 e) de dA ]/cos( i) 

Jo Jo 

This function of the incident angle i has been measured for 
many paints and pigments, while the reflectivity function <j> 
is known for very few surfaces. Since it is difficult to 
relate measurements of I, to measurements of total reflected 
intensity, the device is usually calibrated v/ith the sample 
replaced by a standard of known high reflectivity (e.g. MgO 
or BaSO v powder reflect more than 99% of the incident light 
in the visible spectrum). 
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Figure 8 : Illustration showing quantities appearing in the 
integral for the integrating photometer. 



integral for the Bond albedo. 
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2.1.2.2 PERFECT DIFFUSERS - LAMBERT'S LAW: 

Surfaces made of finely divided powder usually closely 

approximate what has been called a lambertian reflector or 

perfect diffuser. Lambertian emission was first defined for 
black-body radiation and is such that the surface of the body 
appears equally bright from all directions. In the context 
of reflectivity functions we call a surface lambert ian if 

<j> * k cos(i) (the cos(i) accounts for the variation in 

incident radiation). For the highly reflective standards 
mentioned above, we chose k such that all the incident light 
is reflected. 

X'X'lTT */*. 

k J j (1/2 )*s in(2e) de dA - 1 

J O Jo 
k = J/TT 

In addition to the various multiplicative factors shown 
above, a normalized reflectivity function is also used, 
whe re : 


( I, e , g ) = <j> ( i , e , g ) 


yO is called the normal albedo and <|>'(0,0,0) = I. 
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2 .1 .2.3 THE BOND ALBEDO: 

Another integral of the reflectivity function which is used 
is the Bond albedo, defined by astronomers as the ratio of 
total reflected light from a sphere divided by the total 
incident light . 

If the incident intensity is I,, then the ratio of reflected 
light to incident light is: 

r tr/'i 

b(i) 2t rr sin(i) cos ( i) di l/CIiTr 1 ) 

t r/x 

b(i) s in ( 2 i) di 
tt/z. 

<j>(i,e,g) sin(i) sin(2e) de dA di 

2.1.3 THE DISCRIMINANT 7 + 2IEG-(I 3 ’+E" ! ' +G 1 ' ) : 

In this subsection a discriminant is developed which is 
needed in the program implementing the shape-from-shading 
algorithm. This section can be skipped without loss of 

con t in u } t y. 
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The three angles i / e and g, being the sides of a spherical 
triangle/ have to satisfy the following relationships: 

1 + e >/ g y e+g i and g+i e 

It Is often convenient to express these three relationships 

in terms of the cosines 1/ E and G of the three angles. We 

first note that only one of the relationships could fail at a 
time. For example if i+e < g : 

i+2e < g+e i.e. i < g+e and 

2 i + e < g+I i.e. e < g+i 

The angles are all positive and less than tt . Mow assume 
that the condition i + e < g holds, then: 

cos(I+e) > cos(g) 

since cosine is monotonic decreasing for angles between 0 and 
TT. Expanding we get: 

cos(i) cos(e) - cos(g) > sln(?) sin(e) 

Since the right-hand side is positive, the left-hand side 
will be too and we can square the expression. Using I, E and 
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G to stand for the cosines of i, e and g we get: 

( IE -G ) X > <J-I x )U-E l ) i.e. 

I + 2IEG-(I 1 +E 1 +G’') < 0 

We now have to prove the converse i.e. If the angles can 
indeed form a spherical triangle then the discriminant will 

be positive. Since i ^ e + g we have g ^ i-e and similarly 

since e < i+g we have g >/ e-I, so: 

g } |e-i| and similarly 1 ) | g —e | and e > | i -g | 

Applying the cosine as before we get: 

cos(g) - cos(?) cos(e) ^ sin(i) sin(e) 

cos(i) - cos(e) cos(g) ^ sin(e) sin(g) 

cos(e) - cos(g) cos(i) ^ sin(g) sin(i) 

From i+e ^ g etc. we get the same inequalities with the sign 
reversed on the left-hand side. Vie needed to go to the 
trouble of showing that these inequalities hold for absolute 
values of the left-hand side, since it no longer Is 
constraint to be positive. So we have: 



| IE-G| ^ sin(i) sin(e) 
| EG — 11 ^ sIn(o) sIn(g) 
| G I —E | ^ sln(g) sln(i) 
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Multiplying we get: 

| (IE-G) (EG-I) (G I-E ) | ^ C!-I 1 )(;-E I )(r-G 1 ) 

Using one of the two signs for the right-hand side and 
expanding we get: 

(I-IEG)(! + 2IEG-(I 1 +E’-+G a )) >0 

Since |I|, | E| and |G| ^ 1 we have (7 -1E G) 0 and hence: 

1 + 2 IEG-( I^+E^+G 1 ) > y 0 

2.1.4 REFLECTIVITY FUNCTIONS AMD THEIR MEASUREMENT: 

Surfaces where the three parameters i, e and g are not 
sufficient to fully determine the reflectivity are unsuitable 
for this analysis (or at least reduce the possible accuracy). 
Examples are translucent objects and those v/ith n on-is ot rop i c 
surface properties (e.g. the mineral commonly called tiger- 
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eye, hair, thin wax). Perhaps the most important determinant 
of the reflectivity function is the m 1 c ro-st ruct u re of the 
surface (i.e. that structure smaller than the resolution of 
the sensor used in the determination of the reflectivity) and 
different reflectivity functions may apply at different 
magnifications (in addition, at high magnification objects 
become increasingly translucent). It is best then to 

determine the reflectivity function under conditions similar 
to those used in the determination of the shape of the 
obj ect. 

One way to measure the reflectivity function is to employ a 
gonio-photometer fitted with a small flat sample of the 
surface to be invest igated. The device can be set for any 
combination of incident, emittance and phase angles. 

To avoid having to move the source and the sensor into all 
possible positions w.r.t, a flat sample of the surface when 
measuring the reflectivity function, it is convenient to have 
a test-object which presents all possible values of i and e 
for a given g. (The constraints are i+e ^ g, e+g ^ i and 
g + i ^ e ). Use of such an object is greatly simplified by 
using a telephoto lens and a distant source, giving almost 
constant g. It is convenient to tabulate the reflectivity 
versus i and e for each of a series of values of g. A sphere 
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Is perhaps thc« easiest test-object to use if one is willing 
to live with the decreasing accuracy in determining e as one 
approaches the edge. 

One could also have an object of known shape in the same 
scene as the object to be analysed. This solves the problem 
of having to know the source location and the transfer 
properties of the image forming system as well. In some 
cases objects of known shape and surface characteristics 
differing from those of the object under study are useful - 
for example a sphere with specular reflectivity can pinpoint 
the location of the light-sources (e.g. the eyes in a face). 


2.1.5 MATHEMATICAL MODELS OF SURFACES: 

A number of attempts hove been made to predict reflectivity 
functions on a theoroctical basis starting with some assumed 
micro-structure of the surface. White matt surfaces are 
usually finely divided grains of transparent material (e.g. 
snow and crushed glass). White paint usually consists of 
transparent 'pigment' particles (e.g. of S10 2 or Ti0 Z ) of 
high refractive index (1.7 to 2.8) and small size (optimally 
about the wavelength of visible light) suspended in a 
transparent medium of low refractive index (1.3 say). If one 
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Figure 10 : Model of surface-structure. 



Figure 11 : Another model used in the derivation of 
theoretical reflectivity functions. 
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choses a somewhat regular arrangement of suspended particles 
of uniform s i ze and makes some very restrictive assumptions/ 
one can derive a reflectivity function and study its 
dependence on various parameters describing the model of the 
surface. 

Another type of surface is that of a highly reflective 
material (such as a metal) v/here the light rays do not 
penetrate into the material. Choosing a particular type of 
surface depression and a statistical distribution of the size 
of these/ one can again derive a reflectivity function. Only 
a few such models have been studied and little hope exists 
for modelling real surfaces well enough and still deriving a 
closed expression for the reflectivity function. 


2.2 CALCULATION OF IMAGE ILLUMINATION: 

The equations derived in this section are only included for 
reference/ since the program to be described later avoids 
their use by means of a normalization of the image intensity. 
These equations do have their importance however in 
justifying the choice of definition for the reflectivity 
function and in designing optical systems used in 
experimentation with the shape**f rom-shad ing method. 



Page 35 



Figure 12 : Diagram of optical system and quantities needed 
in the calculation of image illumination. 

Let d be the diameter of the pupil of the image forming 
device, f„ its focal length and M the image magnification 
(the ratio of the length of a line in the image to the 
corresponding parallel line on the object). 

Let a portion of the object surface of area dA 0 be Inclined 
at angle e to the line from it to the image-forming device. 

i. 

Its image will have an area of dA* t - dA 0 M /cos(e). 

Let the incident intensity at the object patch be I t per unit 
area perpendicular to the incident ray. Then the emergent 
intensity per unit solid angle will be I t = I, <|>(l,e,g). The 
light captured by the image forming device is l x dA dw/cos(e) 
where dw is the solid angle formed by the cone of angle a. 
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dw = 2TT(J-cos(a/2)) = 4 tt s I n (a/ 4) 

We would like to express this in terms of M and the so-called 
f-number f n : 


fyy * 7/(2 s in (b /2 ) ) 

= /4)+a+V)*( f./d) 1 - (!+!')( fo/d) if ( f o/d) >! 


The f-nurnber usually indicated on the lens-barrel is (f 0 /d) 
or /n/TTTfT/dF . 


ft-U/4) = (I+ M)* ( f„ / d) 1 " 


cos(Q/2 ) 


(7 + M) ( f 0 / d) _ 

y/( Pi V 4) + (7 + M ) x ( f 0 / d )*" 



M *■ 

- 2 TT - 

(4f*+Pi x -7 ) 


if fn> 7 


The intensity per unit area in the image is: 
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l 3 = U [dA 0 /(cos(e) f A )]*[cos(e)/dA 1 dw 
= I x dw/M* - I x 2 rr/ (4 ) if M < 1 

I a - I,^> dw/M a - I, <j> 2rr/(4f*) if M < 1 

It becomes apparent why we chose to define the reflectivity 
function the way we did and also why one might want factors 
of tt and/or 2 in the definition* It should be noted that in 
practice one does not usually employ this equation/ but 
rather normalizes the expressions used. 


2.3 THE IMAGE ILLUMINATION EQUATION: 

This section contains the derivation of the image 
illumination equation and the analytical formulation of the 
shape-from-shading problem. 


2.3.1 PREVIEW OF HOWTO OBTAIN THE PARTIAL DIFFERENTIAL 
EQUATION: 

At a known point on the object we can calculate g. We would 
like to find the gradient (or at least its component in one 
direction) at this point so as to be able to continue the 
solution to a neighboring point. Measurement of the light 
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the geometry of image illumination 
the imaging system. 


on in 
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reflected tells us something about 1 and e. Since only one 
measurement is involved/ we cannot generally hope to 
determine both i and e locally/ but only a relation between 
the two. There are exceptional points where the normal is 
locally fully determined and this is useful in finding 
initial conditions as explained later. This situation is 
contrary to that obtaining in the use of texture gradients 
(see section 2.76) where the gradient Is known locally 
(except for a two-way ambiguity). In obtaining a solution 
from the shading/ only a global approach will work. 

Collapse the two principal planes of the image-forming system 
together/ forming the x-y plane. Let the H-axis coincide 
with the optical axis and extend toward the object. Let f be 
the exit pupil to image plane distance and assume that the 
image and object space refractive indexes are equal. 

Let t be the ratio of image illumination to object luminance 
(can be found from laws of optics - see section 2,2). Let 
a(x/y 7 z) be the incident light intensity (usually constant or 
obeys some inverse square law). Let A(X/y/Z) - t*a(X/y/Z) 

Let r - (X/Y/Z) be a point on the object and r* = (x ,y*,f) 
the corresponding point in the image. 
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b(x',y') is the intensity measured at the image point 

U' ,y '). 


Let I = cos(i)/ E - cos(e) and G = cos(g) . 


We have A ( r) <t>(I,E,G) = b(r') 

rv/ 1 r%/ 


Let p and q be the partial derivatives of 2 w.r.t. x and y. 
We would like to show that this equation involves x, y, z, p 
and q only. 


2.3.2 NOTATION' FOR VECTOR DIFFERENTIATION: 


If A is a vector (3-tuple)/ then A = |A| is the magnitude of 

A 

A. Also let A - A/ A be the corresponding unit vector. 

~ <V IW 

Consider the dot-product A.B as matrix multiplication of the 

^ A* 

1 by 3 matrix A by the 3 by 1 matrix B T (the transpose of B). 

a/ /-> 

Consider partial differentiation w.r.t. a vector as the 3- 
tuple whose components are found by differentiating w.r.t. 
each component in turn. Then for example: 


A 


f\ 


A 
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At times we w i l l also need the partial derivatives of vectors 
w.r.t. vectors. These are defined as 3 by 3 matrices (the 
first row being the result of differentiating w.r.t* the 
first component and so on), then for example: 

f 1 0 o'] 

A « 0 I 0 

\0 o ij 

We will also use partial derivatives of dot-products of unit 
vectors w.r.t. vectors. For example: 


X = A.B and we want X A 
w to 

/v A 

To avoid finding A _ we write A X = A.B and then: 

an* H rv 'v 

AX + A X = A a .B 

A H ~ A 

«v *v /V 


Extending the definition of dot-product in the appropriate 
way we find: 


. [1 0 

A ,B = 0 I 0 B 

~ \0 0 1 J~ 


B 


A X A 



A 

- X A 


X* = (7/A5CB - X A) 
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2.3.3 THE EQUATION IS A FIRST-ORDER NON-LINEAR P.D.E.: 

If the reflectivity function is <J>(I,E,G), the normalized 
incident light intensity at the point r = (x,y,z) is A(r) and 
the intensity at the corresponding image point r' s (x',y',f) 
is b( x *), then : 


A (r) <j>( I,E, G) = b(r^) 

This image Illumination equation is the main equation studied 
here. When finding a solution we assume A(r) and <j>(I,F,G) 
are known and bCr') Is obtained from the image. Vie want to 

/v 

show that the equation is a first-order non-linear partial 
differential equation in two Independent variables, i.e. an 
equation of the form: 

F(x, y, 2 , P, q) = 0 

where p * 2 X and q s 2 y are the partial derivatives of 2 
w.r.t. x and y respectively. Frcm the simple projection 
geomet ry we have: 

x* = (f/z)*r 

/N/ !>J 
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Where f is the exit pupil to image plane distance. We took 
care of image reversal by orienting the x* and y* axes 
appropriately. It remains to show that 1/ E and G are 


functions 

Of X/ V, 2/ 

P 

and 

q. 

An 

i nwa rd 

normal to the 

surface at 

the point r 

is 

n - 

A 

(-P 

/ “Q 

, n. 


Let the 

light-source 

be 

at 

£s 

* (x 

s / y$ / z s 

) . Then t he 

inc i dent 

ray will 

be 

It 

= r 

- Zs 

/ and 

the emergent 


ray -r- = -r. Clearly then: 

_ A A _ /AA A ^ 

I » n . r? , E = n. r„ and G = r;.r« 

/>• *■ a. A. /s- c 

Where the ^ "s denote unit vectors. All the terms thus 

involve only X/ y / 2 , p and q. It follows that we are dealing 
with a first-order non-linear partial differential equation 
in the two unknowns x and y. 

2.3.4 SOME DERIVATIVES NEEDED IN THE SOLUTION: 

When solving the P.D.E. by the method of characteristics we 
will need the following partial derivatives (see section 
2,5)/ which it is convenient to introduce here, following the 
expansion of 1/ E and G in terms of dot-products. Using the 
results developed in subsection 2.3.2 we get: 
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ly. = I r . = (;/r c ) (n - I rO 

A# /V 1 A/ /». U 

In = (I/n )(r - In ) 

/W /s, k O' 

E, = E r t = <I/r t )(n - E r, > 

E„ = O/n Hr, - E n ) 


r * G r + G _ = (1 / r, ) ( r • - G r, ) + (/ / r ; ) ( r, - G 


r• ) 


C 


= 0 


2.3.5 THE EQUIVALENT SET OF ORDINARY DIFFERENTIAL 
E QUAT I ON S : 

The usual method of dealing v/ith a first-order non-linear 
partial differential equation is to solve an equivalent set 
of five ordinary differential equations: 

x = F p , y = F«, , z = pF P + qF<, 

P - “F x - p F^ and q = -F y - qF t 

The dot denotes differentiation vv.r.t. s, a parameter which 
varies v/ith distance along a cha racte r i s t i c strip. The 
subscripts denote partial derivatives. These equations are 
solved along so-called characteristic strips (see [5], page 
24). The characteristic strip are the characteristic curves 
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described earlier (values of x 7 y and z) plus the values of p 
and q on them. 

Since we can multiply the equation F = 0 by any non-zero 

smooth function (x , y, z, p / q) without altering the 

solution surface/ we can obtain a different set of equations: 

x « \ F f / y = \ F, / z * \(p F p + q F, ) 

P = \(**F X - pF* ) and q = \(- F y - q F^ ) 

The solution to this new set of equations wi l l differ only in 
the values of the parameter s at any given point. For 
example if we let 



= ;// Ft ♦ 

F , 1 + 

(PF p 

+ qF^ ) fc ’ 


the parameter 

s g ives 

US 

a rc 

-length along 

the 

cha racte ristics. 

This is 

used 

i n 

the programs to 

be 

described later. 

Of course 

we 

can 

only do this if 

t he 

denominator is not 

zero; at 

s ing u l a r 

points and ambiguity 

edges it will be 

ze ro (i . e . 

Fp 

ir 
ti¬ 
ll 

= 0 and \ oo), 

A 


different choice for X will be used later in the discussion 
of the scanning electron microscope (section 3./). 
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2.3.6 OUTLINE OF PROOF OF EQUIVALENCE OF THE SET OF 
O.D.E.'S TO THE P.D.E. : 

In this subsection the equivalence of the five ordinary 
differential equations to the image illumination equation is 
discussed. The reader who believes the equivalence holds may 
well skip this subsection! 

At a given point (x G , y G , z 0 ) the equation F(x, y, z, p, q) = 
0 represents a relation between p and q. That is, it confines 
the possible solution normals at this point to a one- 
parameter family of directions [53: 

(*p(oO, q(«), -/) 

Increments in the feasible tangent planes thus satisfy: 

dz * p 0 («) dx + q o (ex) dy 
Differentiating w.r.t.cx we get: 

0 = pJCoO dx + qj (cx) dy 

(Dashes are derivatives w.r.t. cx in this subsection). But 
by differentiating the equation F(x, y, z, p, q) = 0 w.r.t.cx 
we also get: 
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FpPoW + F^q'( ) « 0 
Hence: dx/Fp = dy/Fq 

Wha t we need now are similar equations for feasible 
increments in 2 , p and q. First we have: 

dz = p dx + q dy 

in the solution surface (this surface is selected from all 
possible ones by choosing one passing through a given initial 
curve - see later, subsection 2.3.7). Hence: 

d 2 /(pF p + q F s ) * (p dx + q dy)/(pF p + qF<, ) 

= dx/Fp = dy/F<, 

Finally differentiating F(x, y, 2 , p, q) = 0 w.r.t. x and y: 

F x + FjP + F ft ^ + F,q„ = 0 
F y + F, q + F p Py + Fq q y = 0 

Va V* 

but p > " q > = 

Ox ay oyux 

dp/CF^+pF^) = - (p x dx + P v dy)/(p x F p + p y F^ ) = - dx/F p 
dq/CFy+qFji) = - (q y dx + q y dy)/(q x F p + q y F^ ) = - dy/F s 
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dx / F ? = dy/Fc, = dz/(pF p + q F, ) 

= “dp/CF* + pF^) = - dq/(F y + q F^) 

Introducing the parameter s we get the five o.D.E.'s 
mentioned earlier. We have shown that a solution to the 
P.D.E. must also satisfy theses five O.D.E.'s. It is a bit 
more difficult to show that a solution to these O.D.E.'s is 
necessarily a solution of the P.O.E. (see [5], page 28). 
Basically it needs to be shown that the equations for p and q 
produce results which continue to be consistent (i.e. equal 
to the partial derivatives of z w.r.t. x and y) . 


2.3.7 INITIAL CONDITIONS NEEDED: 

To select a particular solution surface amongst all possible 
solution surfaces one needs to specify an initial curve 
through which the solution surface must pass: 

x = x(t), y = y(t) and z = z(t) 


Along this curve we must satisfy: 



z'(t) - p x'(t) +q y' (t ) 

F(x(t), y(t), z(t), p(t), q(t) ) = 0 
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Here the dash represents differentiation w.r.t. t, This pair 
of non “linear equations allows one to find p(t) and q(t) 
along the initial curve (there may be more than one solution, 
in which case there will be more than one solution surface). 
The characteristic strips sprout from this initial curve and 
the solution surface can be described parametrically: 

x = x(s,t), y«y(s,t)/ z « z(s,t) and 
p = p (s , t), q = q (s, t) 

2.4 • SIMPLIFYING CONDITIONS AND UNIFORM ILLUMINATION: 

Since the general equations are fairly complex it is of great 
interest to find simplifying conditions. Some of these are 
presented here, others will be found described in chapter 3. 

I. DISTANT SOURCE: (Collimated source or the object subtends 
a small angle at the source) 
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A „.r: = 0 and for a truly distant source: 

» /v '* 

A r = 0 

A/ 

Replace r- by kr- and let k -> oo then: 

0 , I h unchanged 
0 , E h unchanged 

aw 

(J/r e Hf ; - G £, ) / Gj, = 0 

In addition choosing the 2 -axis along removes further 
te rms. 



2. DISTANT CAMERA: (Telephoto lens or the object subtends a 
small angle at the camera) 

Replace r 0 by kr. and let k-»oo then: 

A# <VV 

I ^ and I * unchanged 

E r = 0 , E n unchanged 

G r = U / r ; ) ( r. - G r { ) , G„ = 0 

In addition choosing the z-axis along r e removes further 


t e rms . 
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3. DISTANT SOURCE AND DISTANT CAMERA: 

I r - 0 , I h uncbanged 
***** ^ 

£ r = 0 , E „ unchanged 
G r = 0 , G h = 0 

Most practical situations are an approximation of this 
case. 

4. SOURCE AT THE CAMERA: 

ri = r e I * E and G = I 

I >. = E h unchanged 
I * * E h unchanged 
G = 0 and G h = 0 

5. DISTANT SOURCE AT DISTANT CAMERA: 

I ► “ E ,. = G r = 0 
I h * E* unchanged, G K = 0 

^ ^ r+ 

Choosing the object to be on the 2 -axis removes further 
terms. This is the simplest possible case. 
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UNIFORM ILLUMINATION: 

Uniform illumination (or an approximation thereof) is 
fairly common and might at first sight appear not to fit 
into our framework. This subsection shows the 
equivalence of uniform illumination to one where a point- 
source is at the camera and a different reflectivity 
f unct i on obta ins , 


The integrals here are analogous to the ones obtained for 
the integrating photometer except that we have constant 
emittance angle rather than constant incident angle. If 
the incident light intensity is I, per unit area oriented 
in any direction, then it Is easy to show that I,/tt falls 
per unit solid angle per unit area perpendicular to it. 
The emitted light per unit solid angle per unit area 
perpendicular to the emitted ray is thus: 

r/t 

I, V^Ce) = 1 1 (7 /tt) [I J <j>(i,e,g)*(J/2)*sin(2i) di dA ] 
/cos(e) 

This is the same situation as if we had a source at the 
camera and a reflectivity function such that: 



<j>(E,E,!> = ^(E) 
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(Except that for uniform Illumination a certain amount of 
se1f-shadowing can occur for non-convex objects) 


2.5 THE FIVE O.D.E.'S FOR THE IMAGE ILLUMINATION EQUATION: 


F(x, y, a, p, q) = A (r) d>(I,E,G) - b(r') = 0 

~ J t-SJ 

We know A(r) and ^>(I,E,G), and obtain b(r^) from the image. 
We need F x , F y/ F ? , F ? and F^. Since £ = (x,y,z) and 

£ * (-p,~q, J) we can get all of these derivatives from F ^ 
and F h . 

•v 

F r = A(r) ^(I,E,G) + A K (r) <|>(I,E,G) - b „ C r') 

Fk - A(r) <|> (I,E,G) 


Let 



(I,E,G ) 

and 


t hen : 




a 


*\ 


Note that a r and a^ are 3 by 3 matrices, the rows of which we 
""" 

computed in a previous subsection (2.3.4). 
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fu/ri )(n - I r ; ) "\ 

'V 1 

(f / r, ) (n - E r,) 

V«*/r t ><jCi- G r t ) ♦ (!/r,)(? t -G r.)J 


f 1/r, -1/n 

J/r, 0 




(J/r, -G/r { ) 


O/rV-G / r e > / 


VW 


Note that this is the product of two 3 by 3 matrices. 
Similarly: 


a 

~ n 


r (l/n)(r ; 
(J/n) (r t 

i 


i 

F. n) 

7 


/-I/n I/n 

-E/n 0 

0 0 



The calculation of b (r') from b ,(r') will be described in 

t* /s/ T *-> 


sect ion 2.7 
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2.6 CAMERA PROJECTION EQUATIONS: 

The projection equations derived here are used in section 
2.7 . 

So far we have assumed the camera to be at the origin 
oriented v/ith its optical axis directed along the z-axis and 
the image-plane x * and y' axes parallel to the x and y axes. 
Moving the camera from the origin introduces only a minor 
change in the equations. If however the camera is oriented 
In a different way, some of the equations become more 
complicated. 

Let R be the orthonormal 3 by 3 matrix which takes the z-axis 
into the optical axis and the x and y axes into the x' and y' 
axes. Then: 


f 

r' = (R r) - 

~ (R r).z 

where z * (0,0,1) is the unit vector along the z-axis in the 
image coordinate system, r' Is the vector from the exit- 
pupil to the image point in the same coordinate system. 

If two images are taken with the camera oriented differently, 
the area recorded in both images wi l l be spatially distorted 
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only. That is, a simple t ran sfo rma 11 on will take the one 
image into the other* 

/ - (R * R ~' £.'> t f ± 

Zx ' (R.RT 1 r' ), z "f, 

Where R, and R z are the two rotation matrices and f, and f z 
the corresponding exit pupil to image plane distances. This 
transformation is useful if one wishes to orient the optical 
axis along r* or r (to simplify the equations for the 

v r 

derivatives) . 

2.7 OBTAINING INTENSITY GRADIENTS: 

To evaluate the derivative F~ (section 2.5) we need b u (r'). 

- X ~ 

b t (r') = b / . r' 

r (w r ’-k/ r~ 

" /V 

f f Rr .2 

r' = Rr-(Rr) --— 

~~ (Rr) .2 ~ (Rr.i)*- 

IV O' « 

f (RrMRi) 

» - (R - - ) 

( Rr) .2 ( Rr) .2 

e- a. 

In the simple case that the camera is oriented properly: 
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R 


Oii) 


f 

f i 0 o) 

/S A 

r z 

r - - 

0 10 


X s /> 

~ r. 2 

A/ 

[o 0 1J 

r . 2 


f 


z 


f* 0 o 
0 10 


f fo 0 x) 
— \o 0 y 

z x \0 0 zj 


f f1 0 -x/z\ 
- I 0 1 -y/z 

z \0 0 0 J 


Written out in full we have: 


(b* , b y , b 2 ) 


( f/ 2 ) [b x /, b y /, - ((x/ z) b x / + (y/z)b y /)] 


b x / and b y , are measured directly from the image. 


Since the intensities measured from the image do not locally 
determine the normal, one might well ask what, roughly, such 
measurements do determine. The components of the gradient of 
the intensity are related to the second derivatives of the 
distance to the surface, while the intensity itself is 
related to the magnitude of the first derivatives. This 
relationship becomes exact for the case of a distant source 
at a distant camera (section 2.5, case 5; see also 3./.2). 

It should be noted that the equation for (section 2,5) 
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also Involves A r , Usually A is fairly constant over the area 

A* 

of the object recorded in the image/ or at least satisfies a 
simple inverse-square equation. 

If A = then A v~ = “2 ( r, 1 /r^ ) r,; , 

* >v ^ 

Where jr ^ is the incident vector/ and r^ is the length of the 
incident vector to the singular point. 


2.8 OBTAINING INITIAL CONDITIONS: 

It would be a great disadvantage if one always required an 
initial curve to start the solution from. Fortunately it is 
usually possible to calculate some initial curve if one makes 
some assumptions about the surface and uses the special 
points where the reflectivity uniquely determines the local 
normal - these points will be called singular points. 


2,8.1 USE OF THE SINGULAR POINTS: 

The singular points are the brightest or the darkest points 
(depending on the reflectivity function). At all other 


points the normal cannot be locally determined. The singular 
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points are points corresponding to values of i and e for 
which the reflectivity is a unique global maximum or minimum. 
These may be either extrema in the calculus sense or at the 
limiting values of the angles. 

This method cannot be used if the surface does not contain a 
surface element oriented in this special direction. The 
points are found by looking for the brightest (or darkest) 
points in the image. 

All we still need to know then is the distance of this point 
from the camera/ but since one is usually only interested in 
relative distances this is not a serious restriction. 
Unfortunately it w i l l be found that the solution will not 
move from these singular points, l.e. x' = y ' =0 . This is 
an indication that the algorithm needs to be informed about 
which way the surface Is curved (convex or concave). To make 
this more concrete assume we have a distant source and can 
thus calculate G at each image point. 
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2.S.2 THE SOLUTION WILL NOT MOVE FROM A SINGULAR POINT: 


Consider the variation of <j> with E first: 

7. If the extremum occurs for 0<E<7 then <j> e = 0 . 

2. If the extremum occurs for E = 7 then n ■ r 0 and hence 
E a = (7/n)(r e - E n> « 0. 

3. If the extremum occurs for E * 0 then n.r„ = 0 and 
E n = (7/n)ftf . That is, xp+yp-z -0 and E p - (7/nr)x and 

^ ^ I 

E<] = (7/nr)y. 

x * <| e (J/nr)x and y = <j> £ (7/nr)y 
2 - <j> E (7 /n r>( px + qy) * <j> E (7/nr)z 

In case 7 and 2 we have cj> E E ? and cj> E E ^ = 0 . 

Now consider the variation of <|> with I: 

7. If the extremum occurs for 0<I<7 then - 0 . 

2. If the extremum occurs for I * 7 then n = r ; and hence 
I K = (7/n)(f. - In) = 0. 

n /vv ^ 
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3. If the extremum occurs for 1=0 then n»r . = 0 and 
I h = (?/n)r:. That is (x-x 0 )p + (y-y 0 )q - (e-2 0 ) « 0 and 
Ip = (!/nr*)(x-x 0 ) and 1^ = (7 /n r^ ) (y-y 0 ). 

x = <j) (7/nr £ ) (x-x 0 ), y * <j> (7/nr 4 -) (y-y 0 ) 

2 * cj> (I/nr. )((x-x 0 )p + (y-y 0 )q) 

= <j) (7 /nr. ) (z- 2 .) 

In case 7 and 2 we have <j> x I p and <j> j 1^ ■ 0. 

Now x = Fp and y = F^ and 2 = pF p +qFc*. 

Fn - A( r) Ad,E,G) 

/V I 

<t> p - <|> T ip + 4 *e e p 

<j>, = <|> T 1 1 + Mi 

So in all combinations of cases 1 and 2 for F and cases 1 and 

• • • 

2 for I we find x = y = 0 and hence also z = 0, therefore: 


That is, the projection of the solution point into the image 
Is not moving as the parameter s is changed. 
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In the case E = 0 we find that even though x and y may be 
non-zero, x' and y* = 0 , 

All that remains is the case I = 0 . Here too x = y = 0, 
when the source is at the camera or if E is undetermined. We 
have found no paints with an extremum for I s 0 where E was 
determined (i.e. the global extremum was not unique). 


2.8.3 GETTING THE INITIAL CURVE FROM A SINGULAR POINT: 

If the surface is convex (or concave) at the singular point 
and we have a guess at the radius of curvature (from the 
overall size of the object for example) we can get around the 
problem of singular points by constructing small spherical 
caps on them. Difficulties will be encountered if this point 
happens to be a saddle point (The presence of a saddle point 
however usually indicates that other singular points exist 
where the surface is either convex or concave). 

Let S be the vector from the camera to the singular point 
(found from its known image coordinates and its distance from 
the camera). R is the estimated radius of curvature and /> 
the distance we decide to step away from the singular point 
(determined in practice by considerations of uncertainty in 
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the position of the singular point and the desired detail in 

A 

the solution)* The known normal at the singular point is N . 

A/ 9 
A 


We construct a spherical cap wi t h center S - RN C . 


Let R* = R'- - r' 

S. = S + <R,-R)N, 

Let X s yxN where y = (0,1,0) 

A -As 

Y = N _x X 

T(t) * p [X cos (2 tt t) + Y sin(2-rrt)] 0<t<1 

a- f a» ~ x 


Points on the initial circle are then given by 


S. ♦ T(t) 

A/ * A- 

Vie also need an initial guess at p and q, so we construct N. , 
(an outv/ard noriral): 

N (t) = R, N + T(t) 

The requirement for an initial guess at the radius of 
curvature is not as restrictive as it might seem, since the 
required accuracy is extremely low. This is because is 
usually very much smaller than R, and hence a change in R 
affects the position of the initial curve very little. Fven 
more importantly the values derived for p and q need not be 
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accurate since they are only used as a first guess in an 

iterative method of finding p and q on the initial curve 

before starting the solution. 

2.9 NON-POINT SOURCES: 

Uniform sources have already been dealt with. Perhaps the 

easiest other case is a circularly symmetric source at a 

distance large compared to the dimensions of the object. 

2.9./ CIRCULARLY SYMMETRIC SOURCES: 

Distant circularly symmetric sources can be replaced by a 
point source after modifying the reflectivity function. One 
merely convolves the reflectivity function with the spread 
function of the source (a bit of spherical trigonometry Is 
involved here). Strictly speaking one should perform the 
same operation with the entrance pupil of the camera since it 
too subtends a finite angle at the object and accepts a 
bundle of light-rays. Since <j> is smooth (except at I M and 
I = /) it will be changed very little except at these points. 
The main change will be that <j> does not tend to 0 as I tends 
to 0, but rather for some negative value of I. Also the 
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specular component will be more smeared out. 



Figure 18 : Illustration of circularly symmetric source and 
quantities used in the convolution. 

^♦ i ■ w— ■mw i ■■ ■■■ ■ ■ ■ i i 

Let the source Intensity be 1(a) per unit solid angle at 
angle a from its center when viewed from the object, 
the new reflectivity function <f*(I,E,G) is: 


tt a»„ 



'in' - a„ 


<j>'(I,E,G)= 11(a) <}>< I",E ,G') a da dv /i 1(a) a da 



0-^0 


t he 
Then 

dv 


Where a D is the total angular diameter of the source. 
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And cos(A) = (cos(g)-cos(1) cos(e))/(sIn(i) sin(e)) 
cos(i') = cos( i) cos(a) + sin(i) sin(a) cos(v) 
sln(^A) = sin(l*) s in (a)/s in (v) 

cos(g') = cos(A+SA) sin(i') sln(e) + cos(l') cos(e) 


2.9,2 MULTIPLE SOURCES: 

When the source distribution is not easily treated as above 
one can introduce a different A^ for each source and replace 
the main equation by: 

^A K (r) <j><I v ,E,G k ) = b(r') 

Difficulties in finding initial conditions will be 
encountered with multiple sources unless they are of special 
kinds (e.g. a point source and a uniform source). 


2.1 0 TYPES OF EDGES: 

Several kinds of edges appear in an image - each with its own 
properties and problems for our algorithm: 
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J. Overlap - (occlusion of one object by another) 

discontinuity in z. The program must detect this or it 
will erroneously continue a solution across such an edge. 

2. Joints - (angular edges on an object) discontinuities In 
the derivatives of z. One cannot continue p and q across 
such an edge. It is possible however to use the position 
of the edge as a new initial curve. This and the previous 
condition can be detected as a step in the intensity 
distribution or from a highlight on the edge. 

3. View edges - special case of 1. , where no joint appears/ 

i.e. the surface is smooth and E tends to 0 as we approach 
it. This is easily detected by the program during the 
calculation of the solution. 

4. Shadow edges - here I tends to 0 as we approach the edge 
and again the program can easily detect this. 

5. Other edge of shadow - if the shadow was bridged this edge 
may serve as a new initial curve. 

6. Ambiguity edges - some are lines of aggregation of 

singular points (on which X ) . The character ist ics 

will not cross an ambiguity edge (see section 3.1.3), 
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2.11 SHADOWS AMD SELF-ILLUMINATION: 

If the single source is not at the camera, shadows will 
appear. Solutions can be carried across shadows since the 
position of the source is known and one can construct a ray 
through the last illuminated point and trace it until it 
meets another illuminated region. Only the coordinates and 
not the local gradient of this new point will be known. It 
is necessary to carry this operation out for all 
characteristics entering the shadow, producing a new initial 
curve at the other edge of the shadow where we can restart 
the solution. In practice care has to be taken because of 
noisyness of the solution. 

Self-illumination is a difficult problem to deal with unless 
the object is convex or its albedo is low (less than 0.4). 
An estimate of the effect of self-illumination can be 
obtained from a consideration of two sem i - i n f i n i te matt 
planes joined at right angles. These ore illuminated from a 
very great distance and such that the incident rays make an 
angle a with one of the planes. Let the reflectivity of the 
surface obey lamberts law and the fraction of the Incident 
light reflected be k. Contrast between two intensities I, and 
I x ‘s usually defined to be: 



SOURCE 
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C = 


I. " I* 
I, + I x 


If we ignore light reflected more than once, we find the 
contrast between the two planes to be: 


C s = tan(a -T X/4) 

While if the self-illumlnation is taken into account we get: 

2 - k 

C x = - tan (a -tt/4) 

2 + k 

Contrast is thus reduced by a factor of (2 - k)/(2 ♦ k) . 
This factor varies from 1/3 to 1 as k varies from 1 to 0. 

Note: the rest of chapter 2 contains some miscellaneous 

items that did not fit In elsewhere. 


2.12 THE INVERSE PROBLEM - GENERATING HALF-TONE IMAGES: 

The inverse problem of producing images of a specified scene 
with shading and shadows is vastly different from the method 
of shape-from-shading. Most programs written for this 


purpose can be used for objects bounded by planes only. The 
main issues of optimization of the calculation of which 
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surfaces are visible to the source and camera respectively 
have been dealt with in some detail in recent work [£] * 
Although the two problems are inverses of one another the 
methods used are quite different. 

An interesting problem of a mathematical nature (and 
incidentally with application to cutting wood-cuts) is that 
of producing curved lines in a plane such that the density of 
lines is proportional to the shading in the image of some 
real or imagined object. Preferrably one would like as small 
a number of 'unneccesary" breaks in the lines as possible, 
i.e. the lines should either close on themselves or leave the 
image. Another restriction one might apply is that the lines 
should not cross (When producing wood-cuts one would most 
likely also reflect some of the surface texture in the choice 
of lines). 

For a special case, a solution is immediately at hand. This 
is the case where we have a distant camera at a distant 
source (section -2.4, case 5; see also 3.1.2) and a 
reflectivity function <j> such that: 

= I = ; / /1 + p l +q l 
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Here the contour tines give a solution/ with no crossing 
tines and no 'unnecce sa ry' breaks. One of the most 
attractive feature of contour maps is perhaps just this fact 
that they provide some shading i n format i on . 


2.73 HUMAN PERFORMANCE WITH MONOCULAR PICTURES: 

Jugding by the popularity of monocular pictures of people and 
other smooth objects, humans are good at interpreting shading 
information. Since they use the same basic information as 
our shape-from-shading algorithm we expect to find similar 
short-comings (see section on facial make-up for example). 
Supposing the human visual system does not use the shading 
information in simple heuristic ways only, one might expect 
that the perception system 'solves' the equations or a much 
simplified form of them. Since this cannot be done locally 
(the way some portions of an edge-finding process might work) 
it is difficult to suggest an elegant and simple mechanism 
and a place to look for it. Presumably it would have to 
involve computational waves travelling outward from the 
s ingular points. 
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2.14 ERRORS AND INCONSISTENCIES: 


It Is difficult to estimate analytically the error in the 
solution because the equations are so non-linear. cj> , b, and 
A cannot be measured to better than 5 or 10% accuracy and 
numerous practical problems such as non-uniform sensitivity 
of the sensor have to be taken care of. 


Only a simple error analysis can be presented here. Suppose 
we wish to determine the effect of varying inclinations on 
how a given error in the input data (intensity in the image) 
relates to errors in the coordinates determined on the 
characteristics. We need to determine the rate of change of 
p w.r.t. b. Consider a particularly simple case, that of a 
distant source at a distant camera (As has been mentioned 
previously and will be demonstrated in section 3.1, the 
equations for this case are particularly simple). Next 
assume that one of the gradient components, q say. Is 0. 


We have 


Then 


b/A 

P 


-_<|>( I) = <|cr/ f I +p*+q x 
= /!/((j>-' (b/A)) 1 - l' 




We need to different I ate p w.r.t. the ratio l = b/A. 
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Pl (4* 1 (i) u ))' 

For both I —* 0 and I —> 1 , the error In p becomes very large 
for a given error In l (since in the first case cj>~' (l) —> 0 
and In the second case <j>" 1 (l) —* J). This is not very 
surprising since in the first case we are looking 
perpendicularly down on the surface and I will vary very 
slowly with p, while in the second case we have near 
tangential incidence and small changes in the angle of 
Incidence (and hence also I) will correspond to large changes 
in p. 

We note that in this rather special case, the error 
contribution to the solution is large in some areas, while 
being small in others where the incident angle is not to 
close to 0 or tt/2. The actual error will also depend on cj)" 1 
and the error in measuring b/A. In a case with less 
restricted lighting conditions the relationship between the 
inclination of the surface and the error-rate will be more 
complex. 

Vie considered the derivative of p w.r.t. I, since It is the 
integral of the error In p which constitutes the error In z 
for any one characteristic. 
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ds 



$ l (s) ds 


Where e(s) is the error in 2 for a given characteristic as a 
function of arc-distance from the singular point, £p(s) is 
the error in p and £l(s) is the error in l. 


In this context one may also want to discuss inconsistencies 
in the solution. If either the lighting conditions or the 
reflectivity function are incorrectly specified, an incorrect 
shape will be calculated. The shape determined may or may 
not violate the requirement of smoothness. If the calculated 
shape is not smooth it can be concluded that the solution (at 
least in some region) is incorrect, and that the given source 
position or the given reflectivity function are incorrect. 
It is easy to give examples of the case where false 
assumption will lead to a smooth solution, as well as those 
where we obtain solutions with discontinuities. 

For simplicity consider a flat, inclined surface (2 = x). 
The characteristics will be straight lines in this plane, 
parallel to the x -2 plane. Modifying the reflectivity of the 
surface to be increasingly darker with increasing x, we 
obtain a new solution which contains characteristics, again 
parallel to the x -2 plane, but curving toward large 2 for 
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large x. This solution is smooth and contains no indication 
of an error. 

If now we apply instead a surface coating which is normal for 
positive y and darker for negative y, we obtain a solution in 
which the inconsistency is apparent. The characteristics in 
the solution for negative y are more inclined than those for 
positive y, and a discontinuity exists at y = 0. 

Using this kind of approach one could determine which kind of 
surface markings are noticeable by an observer (i.e. lead to 
inconsistencies in the solution) and those which merely alter 
the apparent shape. 


2.75 WHAT ARE LIKELY SOURCE DISTRIBUTIONS? 

Since the complexity of the algorithm presented here 
increases with the complexity of the light-source 
distribution and since we only know how to bridge shadows 
cast by one source/ it is important to know which light- 
source distributions occur in practice. First one notes that 
the situations found difficult by humans are almost certainly 
going, to give difficulties to our algorithm. For example/ 
when two sources cast shadows (such as on a road lighted by 
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widely spaced street-lamps) the shape of unfamiliar objects 
becomes difficult to ascertain because of the crossed 
shadows. If the incident intensity varies greatly from one 
image area to another (such as in a lightly wooded forest) 
the tangle of lighted and dark areas makes perception more 
difficult. On the other hand one would expect 'natural' 
conditions to be particularly easy. That is, one point 
source somewhat above the observer (the sun) combined with a 
very diffuse (almost uniform) source (the sky). The diffuse 
source will not throw sharp shadows of its own. The absence 
of either of the two sources makes vision only slightly more 
difficult. 


2.15.1 RELEVANCE TO PHOTOGRAPHY AND GRAPHICS: 

One would expect photographers to have something to 
contribute to this subject and introductory booklets on 
artificial light photography confirm the above conclusions. 
The beginner is advised to use a number of lights with 
different characteristics as follows (Phrases of inexact 
meaning wi l l be placed in quotes): 

J, The main light - The ideal main light is a large spot 
light approximating the effect of the sun. It is usually 



Page SO 


placed 45 degrees above and 45 degrees to the side of the 
subject. Its purpose is to establish the 'form of the 
subject' and fix the ratio of lighted to dark areas. The 
exact ratio is not important but the position of the 
source should result in good shading (which increases as 
the source is moved further from the camera) without too 
much shadow area (in v/hich detail is more difficult to 
pe rceive). 

2. The fill-in light (or axial light) - Its purpose is to 
lighten slightly the shadows cast by the main light and 
approximates the effect of the sky. It is placed near the 
camera to prevent it from casting its own shadows and to 
simulate the effect of uniform lighting (see an earlier 
discussion of uniform i l l um in at I on , section 2,4,6). The 
appearance of shadows within shadows is considered 
extremely 'ugly' and should be avoided since it makes the 
picture more difficult to interpret. The ratio of fill-in 
light Intensity to main light intensity Is usually chosen 
to be about 1 to 3, 

In addition a number of small sources may be used for extra 


e ffeet s : 
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3, The accent light - Its purpose is to enliven the rendering 
by adding highlights and 'sparkle'. It should be a small 
collimated source which can be directed to illuminate 
small sections of the subject. It is placed behind and to 
the side of the subject so that it cannot cast shadows of 
its own. This light can add catchlights (specular 
reflections such as on eyes or metal objects) and bright 
outlines (particularly on hair). 

4. The background light - Its purpose is to 'separate' the 
subject from the background. It illuminates the 
background only, such that the intensity reflected by the 
subject will nowhere match that of the background. This 
ensures that the two can be easily 'separated' - I.e. the 
edge between them will be visible. 

Other hints are that too many lights spoil the effect, having 
the main-light at the camera creates a 'flat' image, shadows 
crossing edges on the subject are to be avoided and that 
light parts of the image draw the attention of the viewer. 
It is interesting to note how much of what is vaguely 
formulated in these introductions to photography can be 
understood from the point of v iew of shading. 
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1 .16 DETERMINING SHAPE FROM TEXTURE GRADIENTS: 

A problem related to that of determining shape using shading 
is that of determining shape from the depth-cue of texture 
gradients. A textured surface will produce an image in 

which the texture is distorted in a way reflecting both the 

direction and and the amount of the inclination of the 

surface. An image of a tilted surface with a random dot- 

pattern for example will be compressed in one direction (the 
average distance between dots is decreased) by an amount 
proportional to the inclination of the surface. Both 

direction and magnitude of the gradient can thus be 

determined - except for a two-way ambiguity. 

In practice it may not always be easy to determine such 
texture gradients reliably because of low resolution of the 
imaging device and scatter/ causing a reduction in contrast. 
Some simple textures may be handled by simple counting or 
distance measurements as suggested above, while more 
complicated textures (e.g. a plastered wall) will need more 
sophisticated techniques, such as two- d i mens i ona l correlation 
(best obtained using the fast-fourier-transform). Some 

experimentation with this technique showed promise, but did 
not supply very reliable gradients and the method was slow. 



Page S3 


The next problem is how to obtain the shape from the texture 
gradients. Starting at some point (whose distance from the 
camera we assume known), we use some external knowledge to 
resolve the two-way ambiguity. We can now take a small step 
in any direction and find the gradient at this new point. 
Continuing in this way we trace out some curve on the surface 
of the object (somewhat analagous to the characteristics in 
the shape-from-shading method, except that here the curve is 
quite arbitrary). 

Let s be the arc-distance along the curve, z 0 the distance to 
the initial point/ and p and q the components of the 
gradient, then: 

z(s) = z 0 «■ 

If one takes small enough steps, one can continue to resolve 
the ambiguity at each step by using the assumption of 
smoothness. This can be done until we meet a point where 
the gradient is zero. To continue past such a point would 
require some external knowledge to again resolve the two-way 
ambiguity. An aggregation of points with zero inclination 

can form an ambiguity edge which cannot be crossed. 
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Clearly we can reach a given point through many paths from 
the initial point. This allows us some error checking, but 
there certainly are better ways of making use of the excess 
information. For that is v;hat v/e have, since we know from 

the solution to the shape-from-shad i ng that only one value is 
required at each point for the determination of the shape, 
while we here have two (the components of the gradient). 
Most commonly when faced with such an excess of information 
on can make use of some l eas t-sq ua re s technique to improve 
the accuracy. Perhaps a relaxation method on a grid would be 
useful (The grid need not be rectangular). 
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3. PRACTICAL APPLICATION: 

3.1 THE SCANNING ELECTRON MICROSCOPE: 

This chapter deals with a few practical applications in which 
the equations simplify considerably. 


3.1.1 DESCRIPTION OF THE SCANNING ELECTRON MICROSCOPE: 

This device uses an electron beam which is focused and 
deflected much like the beam of a cathode ray tube and 
impinges on a specimen in an evacuated chamber [JM. The 
narrow ray penetrates into the specimen for some distance, 
creating secondary electrons along its path (a small number 
of electrons are reflected at the surface). The depth of 
penetration, the spread and the number of secondary electrons 
are all functions of the material of that portion of the 
specimen. The number of secondary electrons which reach the 
vacuum through the surface will depend strongly on the 
inclination of the surface w.r.t. the beam, being least when 
it is perpendicular. 

These relatively slow secondary electrons are then attracted 
by a positively charged grid and impinge on a phospor-coated 
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photomultiplier. In this way a current is generated 

proportional to the number of secondary electrons escaping 
the specimen. There are other modes of operation which do 
not however interest us here. The output is used to modulate 
the intensity of the beam in a cathode ray tube while both 
beams are scanned synchronously in a T.V. like raster. The 
image created exhibits shading and is remarkably easy to 
interpret t opg raph ica l l y. This is quite unlike the normal 
use of optical or transmission electron microscopes which 
portray density and thickness. 

The magnification is easily increased by decreasing the 
deflection in the microscope. The resolution is poor 
compared to the transmission electron microscope because of 
the spread of the beam as it enters the specimen, but the 
depth of field is much better than that of an optical 
microscope because of the very narrow beam (extremely high f- 
number). The higher field gradient on edges causes these to 
be outlined more brightly. This artifact, while appealing to 
people, may be a problem in the implementation of a computer 
algorithm for finding the shape. 

Often the final analysis does not involve exact determination 
of the shape or two stereo-images can be used, but there 
propably are also important cases where the shape must be 
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determined and the stereoscopic method Is not applicable. 
This may be because at the magnification used the specimen 
appears smooth without significant surface detail or because 
it is difficult to line up the second image. Since the 
equations for this case turn out to be so simple it should be 
rewarding to tie a scanning electron microscope directly into 
a small computer. 


3.J.2 EQUATIONS FOR THE SCANNING ELECTRON MICROSCOPE: 

A little thought shows that this is analogous to the case 
where the source is at the camera (or equivalently we have 
uniform illumination); for one thing/ no shadows appear. 
Next we note that at all but the lowest magnifications the 
projection is near-orthogonal. Because of these tv/o effects 
the five O.D.E.'s simplify considerably: 

x = F P / y = F, / 2 = pF p + qF<* 

P s ~ P* “ OF* and q = -F - q 

Now F^=A<j) 1 I n and F r = -b r 
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I - n.z/n = 1/ n (where n = (-- p, - q, J) ) 

I h = (I/n) (5 - In) = (z/n) - (l/n 5 )n 
I ? - U/n 3 )p and I<^ = (f/n 3 )q 

Hence: x = F p = (A <j> x /n 3 )P/ y = F«^ = (A (j» x /n^ )q 

z = (A <j>^/n 3 ) (p x + q 1 ) 
p = -b x and q = -b y 

If <j> T 5 * 0 everywhere, we can change to a new measure s along 
the characteristic by multiplying ail equations by 
X - n 3 /(A ^ 3 ) and we get: 

• • • X 

X = p, y = q , 2 -p + q 

p = b x (n 3 /(A <|) r )) , q * b y (n 3 /(A <|) t )) 

This extremely simple case thus has characteristics which are 
curves of steepest descent (or ascent) , Also note that the 
equation for z does not couple back into the system of 
equations (due to the orthogonal projection) thus increasing 
accuracy. The equations happen to be very similar to the 

eikonel equations for the paths of light-rays in refractive 
media. It may be possible to find ready-made solutions to 
some special cases by using this analogy. 
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We assumed that <f> x j 4 0 ; this is equivalent to assuming that 
an inverse exists which allows us to find I from a 
measurement of the image intensity: 

vp <4>< 1,1,7) ) = I 

Let fCx) = (I- ^(x) )/(2V|/(x) ) 

Then = < I /2 + q*") 

So we can find at each point the magnitude, but not the 
direction of the local gradient. This is very different from 
the method of determining shape from texture gradients 
(section 2.16), where we can locally determine the gradient 
except for a two-way ambiguity. 


3.1.3 AMBIGUITIES AND AMBIGUITY EDGES: 

This is an easy enough example to study ambiguities. 
Consider the two surfaces: 

3 S 

2 = Z + X , H = Z + | X | 

Clearly they cannot be distinguished from monocular v iews 
since their gradient magnitudes are identical: i.e. they 
produce identical intensity distributions in the image. This 



Page 91 


\'t AMBIGUITY 
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Figure 23: A locally determined ambiguity edge. 

o o o 

f = V (X + y -1) 


AMBIGUITY 


bally dete 
/(l+x 2 + (. 


guity edge. 
(l+x 2 + (y+1) 2 ) 
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manifests itself in a slowing down of the characteristics as 
they approach the line x = 0 (alternatively \ DO ). They 
cannot cross this line aggregation of singular points. Mote 
that the characteristics approach this line at right angles 
and that the edge is determined locally, each point on it 
being a singular point. 

A second kind of ambiguity edge can occur parallel to 
characteristics, separating those which can be reached from 
one singular point from those reachable only from another. 
This kind of edge is not locally determined, since a change 
in the surface is possible which removes one of the singular 
points and makes all the character 1stics accessible from the 
other. This can be done without altering an area near two 
given points previously separated by an ambiguity edge. 

Both types of ambiguity edges occur in the general case but 
are not so easily studied there. They divide the image into 
regions within each of which a solution can be obtained. 
Typically most such regions will have one singular point from 
which one may obtain initial conditions (provided one makes a 
decision about whether the surface is concave or convex and 
knows the distance to the singular point). 
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3.2 LUNAR TOPOGRAPHY: 

3.2 .1 INTRODUCTION TO LUNAR TOPOGRAPHY: 

The other very interesting simplification to the general 
shape from shading equations occurs when we introduce the 
special reflectivity function which applies to the material 
in the maria of the moon. This in fact was the first shape 
from shading problem solved both theoretically and in an 
operating algorithm [43. Using the special reflectivity 
function and the fact that the sun is a distant source, it is 
possible (but very tedious) to show that the equations 
simplify so that the base characteristics ( i . e. the 
projection of the characteristics on the image plane) become 
straight lines radiating from the zero-phase point. This 
point corresponds to g = 0 and is directly opposite the sun 
as seen from the camera. Actually this is true only when the 
sun is located at negative z, for positive z (that is in 
front of the camera), the relevant point is the phase 

point, directly in the sun. 
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3.2.2 REFLECTIVITY FUNCTION FOR THE MARIA OF THE MOON: 


The variation of light reflected from the surface of the moon 
with phase and inclination of the surface has been studied 
for a long time. At a given lunar phase g, all possible 
combinations of incident angle i and emittance angle e are 
represented by some portion of the surface. A fairly good 
approximation is the Lommel-SeelIger formula [/]: 

, r. ( i/e) 

Aci,E,G) =--- 

1 (I/E) + X ( G ) 

Where P. is a constant and the function \(G) is defined by a 
table. This formula can also be derived from a simplified 
model of the lunar surface. A slight gain in accuracy is 
possible if H is allowed to vary with G as well. In 
particular Fesenkov [I] finds the more accurate formula: 

r. * <I/E)Cf + cos x («/Z>) 

4<I/E,G) = ----- 

(I/E) + X.(! + tan (<*/2)) 


Where as before: 


tan(o<) = 


G - (I/E) 
/' - G *' 


A recent theoretical model is that of Hapke [33 which 
corresponds fairly closely to the measured reflectivity 
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function. In most of these formulaes we find that for a 
given G, <J) is constant for constant I/E. The lines of 
constant I/E are meridians. 

At full moon, when G = 1 we find that the whole face has 
constant luminosity. This is quite unlike the effect on a 
sphere coated with a typical matt paint where the image 
intensity would vary as: 

yi - (r/R)' L 

Where R is the radius of the image and r the distance from 
the centre of the image. The full moon thus has the same 
appearance as a flat disc if one is used to objects with 
normal matt surfaces. This may explain the flat appearance 
of the full moon. 


3.2.3 DERIVATION OF THE SOLUTION FOR LUNAR TOPOGRAPHY: 

3.2.3./ THE BASE CHARACTERISTICS: 

In the case of pictures taken of the lunar surface from 
nearby (e.g. from orbit) we have the following: 
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1. Distant source (the moon subtends an angle of about ,03 
ml l l i- rad i an s at the sun). 

2. Near point source (the sun subtends an angle of about 10 
milli-radians at the moon). 

3. Camera at the origin. 

4. The reflectivity function is constant for constant I/E. 
This is a property of the material of the maria of the 
moon which has been known for some time. 

We have (using results obtained in subsection 2.3.4): 

U- ~ 0 I n » U/n) (r. - I n) 

E*- = (J/r) (n - E r) E K = (J/n)(r„- E n) 

G r = (11 r) ( r 0 - G r) G * = 0 

~ ~ /V 

Where is a unit vector in the direction from the sun to 

t he moon. 

If I and E depend on some parameter s, while I/E is constant: 

EI S = IE S 
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Since <j>(I,E,G) is constant for constant I/E: 

<j> x l s + (j> g E s = 0 and therefore: 

1 fx + E fe = 0 

If I and E depend on soma parameter k: 

K 1 * + " <k (I * - = <^ r /E)(EI fc - IE k ) 

= 4i E - -«vn 4e 

I j In * 4 E ^ 4 i/ E ^E I K “ IE k ) 

Using some of our previous results we find: 

El- - IE r = - (I/r) (n - E £) 

«L> ^ 

EI h - IE h * (E/n) r 9 - (EI/n)n - (I/n)r + (EI/n)n 
^ /*» ^ 

= U /n) [ ( r.n ) r 0 - (f 0 .?))>] = (J /n ) ( r x r 0 ) > n 
1 

= - ( r x r 0 ) x n 

n v rr 0 

And s ince A s 0 : 

F r “ A " b r 

Vie will ignore F*. for now, mainly because it has an 
looking expansion. 


ug l y 
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F„ = AE^ 


X/G 


nVr, 


( r * r a ) x n 


( r x r 0 ) x n = 

^ v /v 


[x D ( 2-qy )-x ( Zj-qyo )/y 0 (z~px) —y ( z 0 -px 0 ),z(x 0 p+y 0 q)-z (xp + yq)] 


£2 

n 1 - r r, 


». 2E, 

Let X = -AE - where (x 0/ y 0 ,z 0 ) = r 0 . 


Note that A is a constant in this case* 


F„ = X 


F„ = X 


(— 
< Z ° 
Vo 


(7 - 


(7 - 


y 

q-) 
z 

X 

P“) 

z 


x 

-(7 

z 

y 

-(7 

z 


V 

q—) 
z 

X o 

P-) 

z 


Now looking back at the five O.D.E.'s: 


• _ • • • * 

x = F f/ y = Fp / z = pF p + q F,j - px + qy 

p = - F a “ P F^ , q = - F y - q 

Again we can decide to ignore p and q for the tine being, and 
attempt to determine the behavior of the character 1stics, 
Our aim is to show that their projections in the image plane 

are straight lines independent of the scene. The behavior of 

« t « 

y against x is of little help and we next look at the 
projections in the image plane: 
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= x (f/z) and y' * y (f/z) 


(x'/f) = (I/H a )(xz-xs) * (1/z^) [x 2 -x (px+qy )] 


Similarly: 




/ 
z 

X ' 
z 


X X . 

U - -p)x - -qy 
z z 




x e x y X X 

— (J - p~ - q -) - -O - P“ 

z 0 z z z z 


X 

f 


• / 
y 


x r 


( 


x x y 

- -)(1 - p- - q-) 








A 


(--) (1 - p- - q-) 

z 0 z z z 




y 

q- ) 
z 


A 


Now if the surface is not tangent to the ray from the 


E f 0 i.e. r .n f 0 and therefore: 


/V fSJ 


X y 

( J - p- - q-) t 0 
z z 


If in addition \ x / £ t 0, A t 0 and z f 0, then we can 
the two equations : 


yp 

Zo 


_y_ 

z 


• / 

x 


X o 


came ra: 


divide 


o 


H 
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This first-order ordinary differential equation for the base 
characteristics is separable: 

-dy' -dx' 

y e _ _y_l x_o_ _ y/_ 

2 o f H 0 f 

Solving this we obtain: 


y. y x„ x' 

log(-) = log(-) + log(c) 

z» f Ho f 


Let the arbitrary constant c be tan(t): 


7 

sin(t) 




7 

CO s(t ) 


/* 

Xo 

2 o 

V. 





Thus the projections of the characteristics are straight 
lines in the image plane emanating from the point: 



If the sun is behind the plane of the image (z e > 0 - as 

would usually be the case for reasonable illumination and 
avoidance of extraneous light entering the lens) this point 
is called the Hero-phase point, since it corresponds to the 
point in the scene which is directly opposite the sun as seen 
from the camera and hence g = 0, Because of the special 
properties of the reflectivity function of the maria of the 
moon intensity variations in this region are entirely due to 
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non-uniform surface properties rather than shape. It is for 
this reason that this point is not usually included in the 
image but lies somewhat outside it in the image-plane. This 
will prove unfortunate later on when we have to invent 
initial conditions. 

If the sun is in front of the image plane (z 0 < 0), the 
special point Is the TT phase point/ where the image of the 
sun would appear in the image-plane. 


So the obvious simplification to the equations which would 
arise if we let x c s y 0 s 0 cannot be exploited since we do 
not wish to orient the camera in this way. 

We would like to arrange for s, the parameter that varies 
along each characteristic/ to correspond to arc-length. This 
can be achieved by multiplying each of the five O.D.F. s by 
X / whe re : 


(I - p- - q-) 
z z 

Then by choosing constants suitably: 
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XX X, 

- = — = -+ s cos ( t) 

2 f Zo 


y y y„ 

“ = — = — + s sin(t) 

Z f 2 o 


Thus s gives arc length along the characteristics while the 
value of t selects a particular characteristic. 


3 .2 .3.2 THE INTEGRAL FOR z: 


We next turn to 2 which we would like to find without solving 
the messy equations for p and q. 


« ♦ * 

2 = px + py 


x 0 x 


X X [p(-) + q(— - 

2 © 2 2 o 


->1 



Cp cos (t) 


+ q s in (t) 


) 


This is a good place to introduce some abbreviat ions of 
commonly occur ing dot-products : 

*© Vo 

L = (— /—) • (cos(t), sin(t) ) 

2 o 2 o 

x„ y 

M = (p, q) . > 

2© 2 0 
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N = (cos(t), sin(t) ) . (P, q) 


Note that L is predetermined (i.e. independent of the image) 
and that L and M tend to 0 if the camera is pointed directly 
away from the sun (i.e. x 0 ® y 0 s 0) 


x y 

<? - p- - q-> - O-M-sN) 
z z 


z -J 

x = - 

s (J-M-sN) 


and so: 



(J-M-sN > 


We now attempt to express this in terms of measureable and 
calculable quantities (s.a. G, I/E, s and t). Since f 0 

and differentiable it must be monotonic and hence have an 
inverse. That is, given b/A we will be able to calculate I/E 
(G is known at each point). 


r e = (x 0 ,y „ , z „) and r 0 



r = (x,y,z) 


x, 

z (- ♦ s cos ( t ) 

z o 


—— + s s in (t), J ) 

Z o 


Let 
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Then r = 2 O 


n = (-p, -q, 1 ) 
n . r = 2 ( 1 -M-sN) 

ev A/ 

n.r 0 - * 0 (/-M) 


n 




7^ 
+ q 


Let 


r. i. 

T - sL + (—) 


r 0 . r = T z 2 o 

r o • Z 22 o T Zo 

G = ^- = - T = - 

r o r rr 0 Q r g 


So we can calculate G for each point on the cha ract e r i s t i c, 
independent of t and the scene* Next vve attempt to rewrite 
the expression for 2 in terms of I/E: 


n. r. 


I/E = 


2 (1 -M) 


2 Q 


n.r r z(I-M-sN) r 




• Q — 


sN 




J +- 

^ (7-M-sN) 


J 


2 = -“ 


f(I/E) r. A 

-<-; 

Q 2 0 




J 


As mentioned before one can find an inverse ^ to <j) s.t 
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y (b/A ,G ) = I/E 

2 f b r. 1 

z = - - VJ/( — ,G)*(—) - 

s ^ A 2 o Q 

The usual tables for in the case of the maria of the moon 
however are not usually given in terms of I/E and G, but 
rather o< and g. Where : 



tan (o<.) 


G-I/E 



o< is the projection of the emittance angle on the phase 
angle plane. 


2_ T 2 (J-M) 

G “ I/E « -x—- Q —X- 

r Q r D (J-M-sN) 


fl 1. 

r. Q 



sL 


(1 -M) 
(1-M-sN) 



*\ 

J 



f 

(s +L) 

V. 


+ 


N 


(J-M-sN) 



J-G 


2 - 


T 


r 0 Q 


P = sgn ( ?. 0 ) 



De f i n e 
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3.2.3.3 THE INTEGRAL FOR r: 

So far we have been working In the coordinates x', y 
The final result looks neater if we use r instead of 

r = hQ 

r = zQ ♦ z(s+L)/Q 


aN (s+L) 

- Q + 2 - 

O-M-sM) Q 


r 


N 


V 


(J-M-sN) 


Q + (s +L ) 


) 


r = -(rP/Q 1- ) tan(oc) 


an d z. 


Written out more fully / we have: 
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t an (<*) 

r r * X 

s + 2 s L + C—) 

Ho 

The numerator is a fixed quantity for each characteristic, 
the denominator varies along each characteristic (but is 
independent of the scene), white tan(^0 is obtained from the 
measurement of b/A and the known G (using the function 
The given ordinary differential equation for r has the simple 
solution: 



r(s) 

r(0) 


-P 


r 

tan (<*) 



ds 



r(0) is the distance to the point from where the integration 
was started. 


To sum up: as one advances along each characteristic in 
turn, one calculates G, measures b/A and uses to obtain 
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tan(0(), which is then used in the evaluation of the above 
integral. The process is much simpler than the general shape 
from shading algorithm in that the base character i st i cs are 
p re de t e rm ine d straight l ines and only an integral needs to be 
evaluated. It is possible to write the above result in a 
slightly more elegant form, which is the one derived by T. 
Rindfleisch (for h 0 > 0) : 


r (P) 
r(P 0 ) 



tan(oO ds' 


Where s' = fs 


r 0 xr' * -<f/ 2 )(r*r 0 > 



2Ze ( (-~), 

z 2 q 




Z2 0 


) 

2„Z 


y x 

= ~fz e s(-sin(t), cos(t), ——cos(t) --sin(t)) 

Z o ^ 


L = — cos (t ) + 2 

L 


*0^ Vo 

-cos(t) sin(t) +- sin (t) 


z o z 0 


z e> 


N 0 = 




y L x y 

I + —— cos (t)-2-^—* 
2 o ZoZ„ 


cos(t) sin(t)+ —~s in (t) 


= f 2 0 s P 


N ow 
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zxN 0 = -fz a s (cos(t)/ sin(t), 0 ) 

'-s/ /v 

I Z»N, I = sf I z J 

A* 

J N 0 

|ixN 0 | |2 * N 0 | 

*N* ^ A# 

(r'.z) 1 ' = (r.zf = (r.z^/r 1 - = (z/r) 1 '= l/Q 1 - 
The two ways of writing the integral are thus equivalent. 


3.2.4 SOME COMMENTS ON THE INTEGRAL SOLUTION: 

J. The base characteristics are p redete rmined straight lines 
(independent of the image). This makes for high accuracy 
and ease in planning a picture taking mission. 

2. Only a single integral needs to be evaluated, not five 
differential equations. 

3. The primary input is the intensity, not its gradients, 
again making for high accuracy. 

4. Although, as usual, the reflected l ight-intens ity does not 
give a unique normal, it does determine the slope 
component in the direction of the characteristic. 
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J. van Diggelen [23 first noted a special case of this 
when he solved the lunar topography problem for the 
special case of an area near the terminator (line 
separating sunlit from dark areas). The characteristics 
are such that the slope along them can be determined 
locally . The slope at right angles to the 
characteristics cannot be determined locally. 

5. Although T. Rindfleisch [43 did not mention it in his 
paper it is very easy to bridge shadows since each light- 
ray lies in a sun-camera-characteristic plane. Its image 
can thus be traced on the base characteristic until we 
again meet a lighted area. One need not even make special 
provisions for this, but just use tan(*) for grazing 
incidence (intensity = 0 ) in the shaded section. 

3.3 APPLICATION TO OBJECTS BOUNDED RY PLANE SURFACES: 

Since a great deal of image processing these days is applied 
to images of polyhedra one might enquire how one could apply 
this method to such objects. First we note that the main 
features of these objects, the joints (angular edges on an 
object) and edges (where one object occludes another), are a 
stumbling block to the application of our method developed so 
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far. Since we already know that the areas of more or less 
constant reflectivity are plane faces there is little point 
in exploring them. So a completely different approach is 
indicated. Firstly we might check whether a parsing of the 
scene obtained by some other method is correct in the sense 
that the intensities observed for the faces correspond to 
their inclinations (the information on the intensity of the 
faces is usally discarded). It is not clear however what one 
might do if this test fails. Furthermore, the programs which 
reduce the image to a line-drawing and tnen decide which 
faces belong to each object cannot really determine the 
inclinations of the various faces without additional 
assumptions (orthogonality for example). 

One can however find the normals to each face directly from 
the known slopes of the projection of the joints in the image 
and the measured reflectivities. One must know which lines 
in the image are true joints (between two faces belonging to 
the same polyhedron) and which are fortuitous (edges between 
faces of different objects). For each normal we need two 
values. Each intensity gives us one non-linear equation and 
each slope of an image of a joint gives us another. The 
equation from the intensity is of course: 
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A ( r) 


cM I, E, G) 


b(r') 

/V 


0 


Where we know that I, E and G are functions of p and q, There 
will be one such equation for each face. 


Where two faces with normals and intersect, they form a 
joint which will be seen in the image. Suppose two points on 
this image are and JB. Then a vector perpendicular to the 
plane through the joint and the camera is Ax B, We also know 

A/ 

that the joint must be parallel to n * n . But Ax B must be 

t ~ A, 

perpendicular to the joint hence: 


(A x B) , (n y n, ) ® 0 

*** A* f ^ i- 



Figure 25 : Projection of a joint on a polyhedron on 
the image. 
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Each joint contributes one such equation. Next we determine 
how many faces must intersect before we have enough 
information for a solution. Two faces intersecting give us 
two equations from the intensities and one from the slope of 
the image of the joint/ while we need four unknowns. An 
infinity of solutions thus exist. With three faces a 
solution is possible since we have six equations in six 
unknowns. Because of the non-linearity of the eauations 

more than one solution might exist. With a larger number of 
faces we always have at least enough information for a 
solution and at times have more equations than unknowns which 
may remove some of the remaining ambiguities and improve the 
accuracy. In this way too it may be possible to discover 

which joints are really between faces of the same object. 


3.4 FACIAL MAKE-UP: 

When a surface whose photometric properties are taken to be 
uniform is treated so as to change these properties in some 
areas, the apparent shape is changed. This of course is one 
of the uses of make-up. The shape of a face for example can 
be made to conform more closely to what a person thinks is 
currently considered * ideal*. This is achieved by making 
some areas darker (causing them to appear steeper) and others 
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lighter* Areas lightened usually include singular points and 
causo a change in the apparent skin darkness (a noma l i zat ion 

effect) and will change the apparent shape in areas other 
than the singular points. 

These modifications can change the shape perceived v/hen 
viewed under the right lighting conditions. The effect will 
change somewhat with orientation and may at times disappear 
when no reasonable shape would give rise to the shading 
observed. Because of a number of surface oils the skin has 
a specular component in its reflectivity, it is also fairly 
translucent. Both of these effects are sometimes controlled 
with talcum powder. The removal of the specular components 
makes the surface appear more smooth. 
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4 . EXPERIMENTAL RESULTS: 

4.1 A PROGRAM SOLVING THE CHARACTERISTICS SEQUENTIALLY: 

When the solution to the shape from shading problem had been 
found using the inconvenient coordinate system (x', y', z) , a 
program was written which would solve the five O.D.E.'s along 
one characteristic at a time (the equations used are not 
reproduced here). The input data was obtained from the 
image-dissector camera attached (at that time) to the PDP-6 
computer in the Artificial Intelligence Laboratory. This 
camera is a random access device: when given an x and a y 
coordinate it returns a number proportional to the intensity 
at that point in the image. The program first searches for a 
maximum in intensity, constructs a small spherical cap around 
it (to obtain an initial curve) and uses a standard numerical 
method (see subsection 4,7.2) to solve the set of five 
ordinary differential equations. 

The prime data required in this solution is the intensity 
gradient (x' and y * derivatives of the intensity), which is 
obtained from the intensities measured for a small raster of 
points near the current x * and y *, A linear function in x ' 
and y' is fitted to this set of intensities; the 
coefficients of x' and y* are the desired gradients. The 
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size of the raster is chosen to correspond to the step-size 
used in the numerical solution method, so that successive 
rasters almost, but not quite, touch. In this way fair 
accuracy in the determination of the gradients is obtained 
without loss in resolution. 

If the least-squares fit is bad, indicating that surface 
detail is being missed with the stepsize used, or that the 
characteristic is traversing an edge or joint, the solution 
for this c ha rac t e r i s t i c is halted and the solution started 
for the next characteristic. Other reasons for terminating 
the solution are that the characteristic has left the field 
of view of the image-dissector or reached a very dark region, 
most likely a shadow or the background. When either the 
(calculated) incident or emittance angles become very small 
(indicating approach to an edge or shadow edge) or X very 
large (indicating approach to another singular point or an 
ambiguity edge) the solution will also be stopped. 

The data structure here is very simple; just a record of 
various values ( x', y', 2 , intensity, p' and q') for each 
point on each characteristic. The shape so determined can be 
displayed in perspective and stereo on a DEC 340 display. 
The character istics appear as dashed l ines - each dash 
representing a step in the integration (We chose the 
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parameter s so that each dash represents the same arc- 
length). The output can be photographed from the display 
and plotted on a Calcomp plotter. 


4.7 .7 AUXILIARY ROUTINES : 

A number of auxiliary routines needed to be written for this 
program. First the Incompatable Time Sharing System (ITS) 
on the PDP-6 does not support a FORTRAN style arithmetic 
language and all programming was done in assembly language 
(MIDAS). The large amount of arithmetic involved, 

particularly with the inconvenient notation and coordinate 
system used at first, made it imperative to incorporate into 
the assembler the ability to handle arithmetic statements. 

Next we constructed a package of useful routines which 
handles floating point I/O, dynamic array allocation and easy 
generation of display lists for the DEC 340. It also 
includes routines for the standard arithmetic functions 
(SQRT, SIN, LOG etc.) and manipulation of vectors and 
matrices (multiplication, addition, inversion etc.). 
Interrupts, user defined operations and command 
interpret at ion are dealt with as well. Some of the remaining 
routines will be briefly described in the next few sections. 
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4.1.1.1 STEREO PROJECTION AND OBJECT ROTATION: 

Since it is important (pa rt i cu l a r t y during the debugging 
phase) to be able to visualise the shape being calculated/ 
stereoscopic output on the display is provided. 

Let d* be the separation of the eyes, f their distance from 
the display and d 0 the distance from the eyes to the origin 
of the coordinate system (usually chosen to be at the 
singular point from which the solution was started). The 
coordinates of the left-eye and right-eye images of the point 
(x/y/S) are then and (x',y') where: 

x' » (x±dj/2)(f/(E+d 0 )) + d</2 
y' y (f / (z+ d 0 )) 

A pair of lenses is employed to focus on the surface of the 
display while converging on the apparent point (x,y,z). 
Obviously one needs to know the scaling of the display in 
terms of dsiplay units per mm. 

One would like to be able to view the objects from various 
sides and perhaps even have some rotational motion to gain a 
greater perception of depth. To this end the object can be 
rotated around its origin (also offset and expanded in size). 
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Figure 28 : Stereo-projection of an object point. 



Y 


Figure 29 : Definition of Pitch, Yaw and Roll. 
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This was preferred over the more common method of allowing 
the eyes to be moved around in the object space. 


To obtain any orientation with one matrix multiplication, the 
three angles P (pitch), V (yaw) and R (roll) are defined as 
rotations about the x, y and z axes respectively. They are 
applied in that order (order is important because rotations 
a re not commut at ive ), 


A 


f sos ( R) 
* sin(R) 

l 0 


-sin(R) 
cos(R) 
0 



r 

v. 


cos (V) 0 sin(Y) 
0 1 0 
-sin(Y) 0 cos(Y) 


I 0 0 "\ 

0 cos(P) -sin(P ) 

^0 sin(P) cos(P)j 


Using the abbreviation c for cosine and s for sine we have: 


A 


f cR cY 
sR cY 


cR sY sP - sP cP cR sY cP ♦ sR sP"\ 
sR sY sP + cR cP sR sY cP - cR sP 
cY sP cY cP 


The various parameters controlling the object rotation and 
the projective mapping can either be preset or continously 
read in from a number of potentiometers (connected to a 
multiplexor and an A/D convertor) controlled by the viewer. 
While one display list appears, the other Is being calculated 
using the latest set of parameters and v/ill in turn be 
displayed when completed. The parameters are also displayed, 
they are: 
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PITCH, YAW and ROLL (P, Y and R) 

SI EEC or MAG - magnification of the object 

FDIS or DIMG - distance from eye to display (f) 

DOBJ - distance from eye to object (d D ) 

DSEY or EYES - seperation of the eyes (d 5 ) 

For photographic purposes each of the two images In turn can 
be blown up (to account for the reduction In size in the 
camera) and displayed a fixed number of times v/hile the 
shutter is open. Windowing at the edge of the screen is 
automatic and some very simple kinds of hidden line 
elimination are available but net normally used. The same 
stereo display package is used by the later version of the 
program (new SHADE). 


4 . 1 . 1,2 MEASURING THE REFLECTIVITY FUNCTION: 

The reflectivity functions of some paints were measured using 
spheres (large rubber balls) as calibration objects. Both 
camera and source were moved as far away as possible to 
achieve almost constant phase angle g. The image of a convex 
object is especially useful because it contains two points 
for all possible combinations of the incident and emittance 
angles (i and e) for a given phase angle (g). The position 
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of the tight^source is mea s u red, as we It as the distance from 
the front of the sphere to the entrance pupil. The image- 
dissector is focused on the edge of the sphere. 

With the sphere temporarily illuminated from several sources, 
the program finds its exact position and size, as well as the 
difference in horizontal and vertical deflection sensitivity 
of the image-dissector. It is now in a position to calculate 
the points in the image which correspond to given incident 
and emittance angles. For a number of choices of both of 
these angles it then reads the intensity at a small raster of 
points (to reduce noise and the effect of pin-holes in the 
photo-cathode) near these positions and averages them. Since 
there are usually two places in the image with the same 
incident and emittance angle, a check on the data is 
available. The resultant table of values (usually 

normalized w.r.t. the brightest intensity) can be printed and 
the whole process repeated after moving the light-source to a 
new position for a new phase angle. The program accounts for 
such things as change in incident light intensity as the 
light-source gets moved a round. 
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4.1 . 1.3 FINDING THE CALIBRATION SPHERE: 

This subsection and the next deal with details, needed in the 
program for measuring the reflectivity function, which may 
not be of general interest. 

For good accuracy we first need to know the exit pupil to 
image plane distance (the focal length is given). This v/ould 
be easy if one could focus on the front of the sphere. It 
turns out that a simple app rox i mat i on wi l l work in a few 
iterations. At each step one recalculates the estimated 

distance to the edge of the sphere, the estimated exit pupil 
to image plane distance and the estimated radius of the 
sphere, using the previous estimates and the measured radius 
of the image. 

Next we need to find the exact center and radius of the 
sphere from its image coordinates and the known distance to 
its front. First consider horizontal coordinates only. 

x, s f tan(a), x T = f tanCa+b) and x 3 = f tan(a+2b) 

We are given x, and x 3 and wish to calculate Xj, which can be 
done after expanding the tangents. 
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Figure 30 : Determining the exact position of the calibration 
sphere. 



Fjigure 31. Finding points for given incident and emittance 
angles. 
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tan(2b) = (x 3 -x, ) / ( f X -x 3 x, ) 

tan(b) * ( y~] + tan *■( 1 b )' - J)/tan(2b) 

x x /f = Cx, + f tan(b))/(f - x,tan(b)) 

The same formulae are then used to find the vertical position 
y x . Finally we need to find 2 h : 

r * R s tan(b) (y^7 + t a n ^ (b) + tan(b)) 
z = R s f/ /f l t X 1 + y 1 
2 y = 2 (R r «■ r)/R s 


4.7.7 .4 FINDING POINTS FOR GIVEN i AND e: 


Clearly the points for given incident angle lie on a circle 
(on the surface of the sphere). Similarly for points with a 
given emittance angle. These two circles may intersect in 
two, one or no points. One can find this intersection by 
first finding the line along which the planes containing 
these circles intersect. Applying the sine and cosine laws 


we get: 
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Let I = cos(i) as usual and D = |v s | 

D / s in (tt- i)=r/sin(a) b = i - a 

r cos(b) = r(cos(i) cos(a) + sin(i) sin(a)) 

= (r/D) (I /d 1 - - r l (! - I ~) + r (I - I 1 )) 

d = r cos(b) 
v = v, + d v, 

p L '•“'j 

The equation of the plane in which the circle of points with 
given incident angle i lies is: 


v.v< = v .v, = c say (where v - (x,y,z) ) 

One can find a similar equation for the plane in v/hich the 
circle of point with given emittance angle e lies. The 
introduction of an arbitrary third plane allows us to find 
one point v on the intersection of the first two. The line 

A/ 

of intersection of the first two planes must be parallel to 
the cross-product of their normals (let them be v Si and v $z ). 
So the equation of the line we are looking for is: 


(v - v ) = k v, where v = v x v., 

a, ~ 1 ~ t ~S| a- S2- 


The points we are trying to find must also lie on the sphere. 


i . e . : 
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(v - v,) 1 = r 1 

(v + k v. - v ) = r L 

ft. C 

k l v .v + 2 k v. . (v - v c ) + (v — v f ) - r x = 0 

The above equation may have no solution for k, in which case 
no point exists for the given incident and emittance angle. 
Otherwise we can use the two solutions and substitute back to 
obtain the desired coordinates which are then transformed 
into image coordinates. 


4.1.1.5 SOME REFLECTIVITY FUNCTIONS: 

The first paint investigated was a matt white paint 
consisting of particles of S i 0 a. and TiO x suspended in a 
transparent base. Very roughly one finds that the 
reflectivity function behaves like cos ( i) for a given g. 
After playing with polynomial fits for a while, the following 
fairly accurate formula was found by a process of little 
interest here: 


<f(I,E, 


(I+GM2 + G) 


G) = 


r 

i + 


1 + 2 I EG-( I 1 ’ + E 1 '+G 1 ' ) 


16 (I-CO 


Note the appearance of the discriminant discussed in an 
earlier section (2.1.3). The symbolic manipulation program 
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Figure 32: 







Table of 

ref lectivity 

(for 

a 

white 

matt 

paint) 

ve rs us 

I = cos(i) 

> and 

E = 

cos(e) 

for 

G = cos(g) 

= 0.81 

• 

The 


intervals chosen correspond to constant size steps In the 


angles. Note the blank areas for combinations of angles 

which cannot form a spherical triangle (see section 2,1,3) 
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MATHLAB 113] was used to find the derivatives <j) x , and ^ 
needed for the shape-from-shading program. For * reas onah l e/ 
angles the above formula is about 5% accurate, becoming worse 
for extreme angles. The repeatability of this measurement 
was disappointingly low, depending on the depth of the paint 
coat and the conditions of its application. Much of the 
investigation of the behavior of the image-dissector was the 
result of efforts to trace the remaining causes of 
inaccuracy. 

Some other paints and an eggshell showed a matt component 
similar to the above, plus a very strong specular component 
(which is small except near the point for which i - e and i + 
e = g). This component is very sensitive to small changes in 
the surface properties such as can be brought about by 
handling the object. 

The image of a convex object with such a surface will usually 
have two local maxima in intensity. One of these will be 
broad (corresponding to the matt component), the other narrow 
and bright (corresponding to the specular component). These 
may be distinguished by a computer program on the basis of 
just these properties. It would then be possible to start a 
solution from the matt maximum (which is not a global 
maximum) rather than the specular maximum. This might be a 



Page 131 


good idea because of the increased accuracy (for one thing 
the norma t i zat i on of image intensities v/ould be more 
accurate) . 


For the nose-recogn it ion program, a ptaster nose was used 
initially, coated with the matt paint described above. This 
of course was not suitable for the final experiments. The 



A matt paint. 
Lambertian reflec. 
Skin on nose. 


Figure 33 : Comparison of some reflectivity functions. 

restricted lighting conditions described later were chosen 
partly to avoid having to find the full-flegded function o of 
three angles for skin. Since no true sphere covered with 
skin is available, measurements we re taken of the shape of a 
real nose and intensities in an image (of a transparency) 
used to estimate <[>(1,1,7). In this way the non-linearities 

i 

of the photographic process (and they were great) did not 
have to be determined separately. The properties of skin are 
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of course not very uniform and also vary from person to 
person, so no great effort toward accuracy was made. Skin 
has a highly variable specular component, so any 
normalisation had to be done not w.r.t. the brightest point, 
but one nearby. The resultant table of |)(I) versus I lies 
somewhat below the one obtained from the matt paint under the 
same lighting conditions. 


4. 7.7 .6 PROPERTIES OF THE IMAGE-DISSECTOR: 

In an attempt to track down poor results in the first try at 
finding reflectivity functions accurately, the image- 
dissector was investigated in some detail [9], Amongst 
problems found were: 

7. Unequal deflection sensitivity in horizontal and 
vertical directions (differed by 7 1 %). 

2. Twist of image varying with distance from center of 
field of view, 

3. Poor resolution (3 line-pairs/mm - radius of tube SO 


rrm) . 
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Figure 34 : Geometric distortion in image-dissector for a 
triangular raster of points covering the photo¬ 
cathode. (The arrows are exaggerated 3 times.) 
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4. Pinholes In the photo-cathode (about 20 of up to 0.5 
mm in size), 

5. Non-uniform sensitivity of the photo-cathode (varies 
more than 30%). 

6. Fairly long settling time of the deflection coils. 

7. A large amount of scatter, which reduces the 
contrast by almost one-half and causes intensities 
measured on the image of a uniform square on a dark 
background to vary by 2 0%, depending on how close to 
the edge the measurement is taken. 

Some of these difficulties are inherent in the state-of-the- 
art of these devices, others were repaired. In any case, it 
was possible now to think about how to improve the program to 
be more insensitive to these shortcomings. 

The program for finding reflectivity functions using spheres 
as calibration objects was sensitive to the (at that time) 
severe deflection inaccuracies, since the emittance angle 
varies rapidly near the edge of the sphere (this effect could 
be reduced with a parabolic t es t-ob j ect) , A calibration 
table was created by another program in which are recorded 
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the image-dissector coordinates of a rectangular raster of 
equally spaced points on the photo-cathode. Also recorded is 
the sensitivity of the photo-cathode at each point. A simple 
interpolation routine can then be applied to coordinates sent 
to the image-dissector to counteract the distortion and, 
similarly, the intensity values returned can be corrected. A 
more convenient triangular raster of points covering the 
whole photo-cathode was later established. Adjustments to 
the image-dissector eventually reduced the distortion by a 
significant factor and use of this table was no longer vital, 
although it did improve accuracy. 


4.7.2 NUMERICAL METHODS FOR SOLVING THE O.D.E.'S: 

The five O.D.E's were at first solved using a well known 
Runge-Kutta method [7, page 2723. The idea is that at a 
given point we can calculate the derivatives of the five 
variables (x', y', z, p' and q') w.r.t. the parameter s. 
Using these we take a half-step forward (increment s by h/2 
and calculate new values for x', y', z, p* and q' as though 
derivatives higher than the first where zero). We then 
calculate the derivatives at this new point and take the same 
step (now using the new derivatives, which will differ 
slightly from the previous ones). We next take a full step 



SHADE SCRIPT 70:05:18 08:15:53 P*ze 13i 

PITCH= 0.400 VAU= 0.000 ROLL= 0.000 SIZEC= 2.0 
FDIS= 118.7 DOBJ= G08.2 DSEV= 5S.0 




Stereo pair of solution produced by old SHADE. 




Figure 36 : Same solution viewed after a slight rotation 
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(increment s by h) The final full step is taken using a 
weighted average of the four derivatives found In this way. 
Written out in symbols this becomes: 

Let h be the step-size (for the parameter s) 

And V = (x',y' /2 ,p',q') 

Also let the equations for the derivatives be: 

Y' = F(s,Y) 

rv n 

(In our case F is actually Independent of s) 

Denote Y(s B ) by Xu then at the n ^ step: 

b - h F(s ,' x»> 

= h F(s h +h/2, Y h + K i / Z) 

Kj = h F(s„+h/2, Y„+K,./2) 

K v = h F(s^h / Y h+ jS 3 ) 

J* + , = Xh + (//6)(K, ♦ 2K 1 + 2K 3 +K V ) 

This method is easy to start (requires no previous values of 
Y) and stable, but requires four time-consuming evaluation of 
the derivatives per step. For this reason various predictor- 
modifier-corrector methods [7, page 1941 were tried and the 
simplest was found to give adequate accuracy: 
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£ *+< 

* I* 

♦ 

Z h F(s„, Y h ) 

Sn + . 


- 

(4/5 ) (P M - C„) 

£am 


+ 

(h/2 ) ( F(s k , M h+1 ) ♦ F(s h , Y„>) 

£h +1 

£ *+1 

+ 

U/5)CP„ + , - C h „ ) 


P, M and C are the predictor/ modifier and corrector 

(V 

respectively. This method is stable and requires only two 
derivative evaluations per step, but is not se l f-start i ng. 
The Runge-Kutta method was retained for the first step in the 
integration. Stability and accuracy were not serious 
concerns since the noise in the data input contributes far 
more to errors in the solution. 


4.J.3 ACCURACY OBTAINABLE: 

Under optimal conditions (using the methods described to 
cancel out most of the distortion and non-un i form i ty of 
photo-cathode sensitivity) the program was allowed to scan a 
sphere of 100 rrm radius, A sphere was then fitted by an 
iterative least-square method to the data points found. The 
data points nowhere deviated from the fitted sphere by more 
than 10 mm, and by less than 5 rrm except near the very edge 
of the image. Such accuracy will not usually be obtained 
because of non-uniformity in the paint/ shortcomings of the 
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sensing device etc. For many purposes however less accuracy 
is quite acceptable and for object recognition in particular 
a more important criterion is most probably that similar 
objects are distorted in similar ways. 

4.1.4 PROBLEMS WITH THE SEQUENTIAL APPROACH: 

It soon become apparent that solving the characteristics 
sequentially had many disadvantages in the general case, even 
though it works well for lunar topography. The first reason 
is that as the characteristics spread out from the singular 
point/ they begin to separate and leave large portions of the 
image unexplored/ obtaining only a very uneven sampling of 
the surface of the object (This is no problem for lunar 
topography since here the solution is not started from the 
singular point/ but at a place where the spread of the 
characteristics is small). 

With a more parallel approach new characteristics can be 
interpolated as we go along (and some deleted as they 
approach too closely). 

Next we find that the base characteristics (projections of 
the characteristics onto the image) may sometimes cross. 
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IMAGE AREA 




Figure 37 : Comparison of spread of characteristics for 
typical solutions in case of lunar topography 
and general case. 

• —.. 

This is not possible if the solution was exact, since it 
indicates that the surface is double-valued or at least that 
its gradient is double-valued. Characteristics may converge 
or diverge from a (singular) point however. Crossing of 
characteristics is really symptomatic of another problem 
which was touched upon when proving the equivalence of the 
five O.D.F.'s to the P.D.F. : The differential equations for 
p and q must continue to give consistent results with the 
surface calculated - this does happen if the solution is 
exact, but cannot be hoped for with the noisy data obtained 
from the image. What one would like to do is continuously 
monitor whether the current p and q match with the slopes 
obtained by first differences from points on the current and 
neighboring characteristics. This is not possible if the 
characteristics are solved separately and are spreading apart 
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a s V;e ^ • A later section {4*2,2.]) will explain a method 
used to cont inously adjust p and q to rerin in consistent 
(derived from the method explained earlier for finding p and 
q on the initial curve). 

At a very minimum, to avoid embarassrnent one would like to 
detect when two characteristics approach one another and stop 
one before they cross. This is easy if the solutions are 
carried along in parallel, but involves lengthy comparison 
tests otherwise. 


4.2 A PROGRAM SOLVING THE CHARACTERISTICS IN PARALLEL: 

Once it had been demonstrated that the equations wore correct 
and a numerical solution possible it was decided to write a 
second program which would sample the surface of the object 
more evenly by interpolating new characteristics when needed. 
Less attention was paid to accuracy in the solution while 
attempting to be less sensitive to various noise-effects. At 
the same time an effort to. find a more convenient coordinate 
system produced the much shorter notation and resultant 
equations described in chapter 2. The solution is achieved 
by taking all characteristics one step forward at the sorxB 


t i me . 
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4,2.1 THE BASIC DATA STRUCTURE: 

The values stored for each point (x.. y, 2 , intensity, p, q 
and pointers to the previous point on the same 
characteristic) are here arranged not by characteristic but 
by 'ring*. A ring is a curve of constant arc-distance from 
the singular point - i.e. the n points on all the 
characteristics form one ring (arranged in countet clockwise 
order of the corresponding image points). The complete data- 
structure is made up of a number of rings, the first of which 
is the initial curve. As before, individual characteristics 
may stop for a variety of reasons (s.a, crossing an angular 
edge) and this causes breaks to appear in the current ring. 
The break is indicated by a point having a negative 
intensity, the value being a code for the cause. Seme rings 
thus represent closed curves (e.g, the initial curve) and 
others more distant from the singular point are broken into 
sections, the final ring having no active point on it (i.e. 
positive intensity). Scavenger routines are usually invoked 
at each solution step and amongst other tasks, compress 
series of dead points (i.e. negative intensity) into one, 
since only one is needed to mark a break in the ring. 

As we have seen one of the main inducements for using the 
parallel solution method is to allow interpolation of new 
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characteristics - this is one of the reasons why the number 
of points in a ring may change from one to the next and why 
each point has to have a pointer into the previous ring, 
indicating which element is its predecessor in the same 
characteristic. This pointer is ~1 if no previous point 
exists (e.g. on the initial curve or the first point in an 
interpolated characteristic). 

We have seen how character ist ics may be terminated causing a 
break in the ring; it is also possible for a characteristic 
to disappear/ without causing a break/ when two 
characteristics approach too closely. In addition a break 
can reclose if the points on either side of the break are 
within the step-sise (and pass the crossing-test explained 
later). With all of this In mind it becomes clear that the 
dat a ~s t r uct u re can at times look pretty confused and this has 
to be remembered when defining a function which interrogates 
the neighbors of a point (s.a. some sort of difference 
appp roximation). 

It was decided to use as data only the coordinates and the 
slo pe at each point/ because this was sufficient for t he uses 
to be made of the data and also was easily available. For 
some uses more complicated surface descriptors may be in 
place, such as the rational function approximations for each 
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surface-section described by Coon [1 0]. Usually the 
increased complexity Imposed by such an approach can be side¬ 
stepped by rather using a smaller stop-size to obtain a finer 
grid. 

It should be noted that the user of constant size steps along 
the cha racte r is t ics may produce difficulties on complex 
objects. For even with smooth surfaces the curves of 
constant arc-distance from the singular point (the rings) may 
have cusps. This invalidates the use of difference methods 
on points along these curves (s.a. are used in subsection 
4,2.2.1 and 4.2.2, 3). No difficulty was experienced with 
images of the objects we experimented with. An alternative, 
which would circumvent this problem, would be the use of 
steps traversing a constant increment in intensity. This 
would turn the rings into contours of constant intensity. 


4.2.2 EXTRA PROCESSING POSSIBLE: 

4.2.2.1 SHARPENING - UPDATING P AMD q: 

We have already described how one can obtain p(t) and q(t) on 
the initial curve by solving the set of non-linear equations 
(subsection 2.3.7): 
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P(t) X fc (t) + q (t) y t (t) - J^Ct) = 0 
A( r) A( I,E,G ) - b( r' ) = 0 

** I /N J 

In the proof that the solution of the five ordinary 
differential equations is also a solution of the original 
partial differential equation, it was stated that the two 
equations for p and q do in fact continue to give the 
derivatives of z w.r.t, x and y. When solving a difference 
equation approximation from noisy data we can expect the 
solution for p and q to become progressively more inaccurate. 
Yet the above pair of equations must hold on any path along 
the surface of the object. In particular one can use them on 
the curve defined by one ring to determine values of p and q. 

For the initial curve we had the additional difficulty that 
the two equations might have more than one solution and we 
selected one on the basis of some external knowledge Ce.g. 
that the object is convex near the singular point). We have 
already assumed that the object is smooth and therefore we 
will have fairly good values for p and q and cannot get into 
this difficulty at non-singular points. Even a simple 
Newton-Raphsen method will suffice to get us more accurate 
values of p and q. 
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Let g(p,q) - p x* + q y*. - z t 

h(p,q) = ^(I,E,G) - b/A 

And suppose: p;(p + £p,, q + £q) = h(p+ip, q+^q) = 0 


Then ignoring other than first-order terms we have: 





q ; 


r g (P/q) 
^h (p / q) 


That is : 



He re x ^ an d y t 

have 

to 

be 

es timated 

f rom 

d i f fe rence 

app roximations . 

One 

may 

not 

want to 

app l y 

the full 

correct ion (£p/ £q) 

. More 

than 

one iteration 

will not be 


required since p and q are very close to the correct values* 


4.2.2.2 INTERPOLATION AND CROSSING TESTS: 

When the separation between two neighboring points In a ring 
becomes greater than 1.5 times the step size along the 
characteristic/ a new characteristic is interpolated. Its 
X/ y / z, p and q values are set to the average of its 
neighbors while the backward pointer is set to -I. A more 
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complicated interpolation method can also be used which 
constructs the line of intersection of the tangent planes at 
the two neighboring points. It then finds the point on this 
line closest to the two neighbors and finally uses a point 
half-way between the point determined previously by the 
simpler method and this new point (This, for small angles 
between the tangent planes, is accurate for a spherical 
surface). This method does not however add significantly to 
the accuracy of the solution. 

If two neighboring points in a section of a ring come closer 
than 0,7 times the step-size, one is deleted (It is important 
that this factor be less than 0,75, that is, one half of the 
factor used in the interpolation decision, or succesive rings 
on a flat region w ill have points interpolated on one step, 
only to be removed on the next, with consequent loss of 
accuracy), 

Finally one wants to stop neighboring characteristics from 
crossing over each other. Consider the two points a and b on 
one ring and their successors c and d on the next. The test 
consists of checking whether c is to the left of the directed 
line through bd and whether d is to the right of the di rected 
line through ac (Both tests are needed). If either fails, 
the corresponding characteristic is terminated, causing a 


P a^e 1U 


d 



Figure 38 : The four points used in the crossing test. 



Figure 39 : The five neighbors used in determining the 
intensity gradient at P. 



Figure 40 : Covering the image with the rasters of points 
read for each solution point. 
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break to appear In the ring at that point. The test is 
equivalent to checking whether the line segment cd falls in 
front of the line segment ab (and does not cross it). This 
test is applied across short breaks in rings as well, to stop 
neighboring section of the ring from crossing over each 
other. 

Care has to be taken if the sections of a ring left all fall 
on one side of the singular point / since the break then 
actually encompasses an arc of more than TT and crossing tests 
applied across it will invariable terminate more 
characteristics on either side of it. This can be avoided if 
the crossing test is not applied to points whose images fall 
too far apart (in terms of the projection of the current 
step-size ), 


4.2.2.3 OBTAINING GOOD INTENSITY GRADIENTS: 

To be more noise-immune than the previous program, a better 
way had to be found to obtain intensity gradients. Rather 
than use the intensities at a small raster of points to 
estimate the local gradient, it was decided to use a 
difference approxImation from intensities measured at 
neighboring points. Using as many as possible of the 
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intensities of the point itself and its five i nm e d i a t e 
neighbors, we can apply a simple leas t-squa res method to 
estimate the gradient. Some of the points may not exist as 
explained previously and the characteristic is terminated if 
less than three points are available or only three which are 
nearly colinear. Suppose the coordinates of the points are 
( x k/Y^) (image coordinate system) and the intensities are b K . 
We wish to find b 0 , b x ' and b y >, to minimise the following 
exp ression: 


♦ S'yf + b o - v 

k 


This happens when: 




* * 


x^y K 


x kVk 

* x 

y t < 


V 







b K>k 


\-> b K 


J 


From b x / and by, v/e can find b x , b y and b by using the 
camera projection equations of an earlier section (2.7). 


For good noise-immunity and some ability to detect surface 
detail indicating that the solution is invalid, the intensity 
for each solution point is not read from only one image 
point. Small tilted rectangular rasters of points are 
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established around each point of the solution. The one axis 
of the rectangle is parallel to the base characteristic at 
that point/ and the s i 20 is adjusted to correspond to the 
projection on the image of a square on the object of side- 
length equal to the step-size. The intensity recorded for a 
solution point is the average of the intensities read for the 
points in this raster and the r.m.s»/average is used to make 
the edge-crossing decision. The rasters of all the points in 
the data-structure almost/ but not quite/ touch and taken 
together almost cover the total area of the image explored. 
This insures that the data is not much affected by pin-holes 
in the photo-cathode of the image-dissector and that edge 
crossing can easily be detected/ without reducing the 
resolution . 

Both this program and the one discussed in section 4.1 spend 
more than half their time accessing the image-dis sect or. 
Between 20 and 100 intensities are read for each point In the 
solution/ and each access takes about .2 to 1.0 milli¬ 
seconds. A complete solution requires from 1 to 5 minutes 


of real time. 
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4.2.5 A DOZEN REASONS TO TERMINATE A CHARACTERISTIC: 

This is a good place to summarise the reasons for terminating 
the characteristics. The values printed near the end of a 
cha racte r ist i c (derived from the negative intensity code 
discussed earlier) can be used to index this table. 

1. The cha racte r is t i c has moved out of the field of view 
of the image-dissector, 

2. The r.m.s./average for the intensities read in the 
raster has become too great, indicating overlap of two 
objects or an angular joint on one object or some 
surface detail that is being missed. 

3. The intensity has become too low, indicating a shadow 
region. 

4. X is too large, indicating approach to either another 
singular point or an ambiguity edge. 

5. There are too few neighbors to construct a good 
estimate of the local intensity gardient. 
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6. is too snail. CX is the Jacobian of the image 
transformation from 2 X and z y to 2 *' and 2 y '. This 
transform becomes singular when CX = 0, In most cases E 
will become too small before this happens. 

7. A new point was interpolated but both its neighbors 
were terminated before it could get anywhere. 

S. I too small - indicating approach to a shadow edge. 

9. E too small - indicating approach to an edge of the 
object. 

10, This c ha ract e r is t ic crossed over a neighboring one. 

11, It was discovered that this point has a backward 
pointer to a nonactive point. This is really an error 
condition and shouldn't normally happen. 

12, The intensity is equal to or greater than that measured 
at the singular point, indicating another singular 
point or ambiguity edge. 

Note that several of these conditions are redundant to ensure 
that even with an inexact solution at least one will fail at 
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the right place. 


4.2.4 OPERATION OF THE PROGRAM: 

4.2.4./ THE INTEGRATION PROCESS 

First the program needs to be given such parameters as the 
position of the light-source, the distance to the object/ 
focal length of the lens and the step-size to be used in the 
integration. It then proceeds to find a point of maximum 
intensity (for some reflectivity functions one needs to 
search for a minimum). This search can be directed to allow 
a choice of one of several possible maxima. The program then 
assumes that this point of maximum intensity is a singular 
point and that the object is convex at this point (in some 
cases we would like to assume it to be concave). After 
constructing an initial curve (a small circle) around the 
singular point, it proceeds to read the intensities at the 
corresponding image points. The non-linear equations for p 
and q on this curve are then solved iteratively. 

All intensities are normalized w.r.t, the intensity at the 
singular point unless the surface has a specular component. 
In the latter case, the intensities on the initial curve are 





FIGURE 41C 































































































Figures 42 A, B, C : Stereo-pairs of solutions produced by 
new SHADE for disc-shaped, spherical and bullet-shaped 
objects (actually spheres with make-up applied). 
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Figures 43 A, B, C : Stereo-pairs 
previous figures, rotated 90°. 



of same solutions as in 
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used to establish a normalisation value (The specular 
reflectivity is too variable for use in norma l i sat ion) . It 
is assumed that the initial curve has been chosen large 
enough fo fall outside the region of strong SDecular 
reflection. 

For each step in the parameter s, the following procedure is 
t hen ca r r ie d out : 


1, For each point calculate the normal (n) / the incident 
vector ( r c - ) and the emittance vector ( r a ) . From these 

^ -v- *■ 

obtain the derivatives l n , E K and G h . 


2. Calculate i , |) G and hence 


* * 


3. Then obtain F , F^ andX. 


4. Add (^x / £y, &s) to (x / y / s) to get the point on the next 
ring for each cha racte r 1st ic. Here (£x, i y, £z) = 

\(F f/ P Fp +q F^ ) « 


5. Interpolate new points where the points in the new ring 
are too far apart and delete points where they are too 
close together. Produce breaks where cha ract e r i s t i cs 
have crossed over adjacent characteristics. 
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6 . Now read the intensities for all the points. Terminate 
those characteristics with points of very low intensity 
or high r.m.s/average. 

7. Calculate by, by for all those points for which enough 
neighbors exist. From these values obtain b*, b y and b^ 
by the projection equations. 

8. Now use n, r- and r„ to calculate I P , E y and G^.. 

ft* ft. ^ ft. V ft. 

9. Next use <j> x , <j> e and ^ to calculate <|y. 

10, Then obtain F x , F y and F^. 

11, Add (£p, Sa) to (p,q) to obtain p and q for the 

un inte rpolated points on the new ring. Here (£p,£q) = 

\( (-F x -pF^ ), (~ F y - q F^ ) ) . 

12, Interpolate p and q for the new points, 

13, Sharpen up the values for p and q on all points in the 
new ring. 

14, Garbage-collect various items, such as series of points 
with negative intensity. 
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15, Stop if no points with positive intensity remain* 

It should bo apparent where the various tests for terminating 
the characteristics fit into the above schema* The simple 
Euler method for solving the differential equations could be 
replaced by a Runge-Kutta method with increases in running 
time of a factor of two, but little improvement in accuracy. 
The sharpening method, on the other hand, is very cheap and 
contributes substantially to accuracy. 


4,2,4,2 OTHER PROCESSING AVAILABLE 

As explained before, the data-st ruct ure is displayed as it is 
generated and can also be viewed from different angles when 
completed. In addition a mode exists where the mapping from 
three-space to the display surface is not performed by the 
projection explained earlier, but a simple map from a 
rectangular area on the image- d i s sect or to a rectangular area 
on the display surface. This is particularly valuable for 
overlaying the solution on an intensity modulated display of 
wha t appears in the image. This aids greatly in deb ug g i n g 
since it is easy to pinpoint such problems as starting the 
solution from an inappropriate maximurn in intensity. 
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A number of other displays can be produced to aid in setting 
up the image-dissector. Prodigous amounts of detailed print¬ 
out can be generated during a solution process and a more 
parsimonious listing of the final data is available. It is 
possible to substitute synthetic data (with selectable 
amounts of noise) for the image-dissect or input as a 
repeatable way of checking out the program and to tide over 
those days when the image-dissect or is being repaired! The 
data can be written to and read from the disk and tape.' 

The stereoscopic display has to be viewed with an appropriate 
pair of lenses which are not always handy. For this reason a 
routine was provided which produces a contour map from the 
data. This map is produced by first listing the 

intersections of all the lines in the data structure (from 
point to point in each characteristic/ as well as from point 
to point in each ring) with the selected contour planes. The 
intersections are then sorted on contour plane. ’Within each 
contour plane the following process is applied repeatedly 
until no points are left: 

Pick a point and find the closest neighbor within a 
'reasonable' distance. 'Reasonable' distance is defined 
to be 1,5 times the step-si?.e used in the solution. Now 
another point is selected closest to the new point also 
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within a reasonable distance and so on until no more can 
be found. The point chosen at each step may not be the 
first in the chain so constructed (which would close the 
loop) unless no other points are available. Also the 
line-segments connecting sucessive points may not make 
angles of more than jT/2 with one another. The points 
are removed from the data as they are used in generating 
the contour except the very first point (to allow for 
the eventuality of closing the contour). 

The distances are usually weighted with the dot-product of 
the new segment with the previous segment/ to give preference 
to continuation of contours in the direction of the last 
segment used. The method generates good contours where the 
data is complete and smooth, and does fairly well otherwise. 


4.2.5 INSENSITIVITY TO IMPERFECTIONS IN THE SENSOR: 

This program is not quite as accurate as the one that solves 
the characteristics serially (mostly because of the simple 
method for solving the differential equations numerleatly), 
but vastly superior in its behavior v/hen faced with noisy 
data. Most of the improvement is due to the better way of 
obtaining intensity gradients and to the use of the lateral 
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FIGURE tfA 
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FIGURES 46 A,B,C 
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SHADE 314 70:05:07 13:25:52 FUM 2 

PITCH= 0.000 YAW=-1.300 ROLL= 0.000 MAG= 3.0 
DIMG= 114.0 DOBJ=1000.0 EYES= 68.0 



for the plaster object. 
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connection between the characteristics* The difference 
approximation for the intensity gradients uses a support area 
about six times as large as the one used by the least squares 
approximation of the first program. 

Distortions in the imaging device 'merely' produce 
distortions in x and y, while non-uniformities in the 
sensitivity will affect p and q and hence 2 . The only effect 
of low resolution will be that some edges will not be noticed 
and the solution erroneously continued across them. 


4.3 A NOSE-RECOGNITION PROGRAM: 

To illustrate one use of the shape-f rom-shad ing method / it 
was applied to a simple recognition task. Although there is 
great interest in f ace-recogn i t i on [12] (partly because there 
is a practical use for it), it was decided to tackle a sub¬ 
problem - that of nose-recognition. In principle, face- 
recognition could be carried out by repeating the process 
explained here for not only the nose, but the chin, forehead 
and the tv/o cheeks. T ran spa renc i e s of noses, rather than 
real noses were used because they are always ready and do not 
move during the minute or so it takes to determine the shape. 
To avoid having to determine the reflectivity function for 
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skin as a function of alt three angles, special lighting 
conditions were employed. The light-source was placed near 
the camera and the reflectivity function as a function of the 
incident angle determined from the tran sparencies taken. 
This meant that no separate determination of the non- 
linearities in the photographic process was needed. 


4.3.7 MODIFICATIONS TO THE BASIC PROGRAM REQUIRED: 

A few minor changes and additions had to be installed in the 
main program for this task. Most prominent amongst these is 
the procedure used to normalize the intensities read from the 
image. Because of the strong specular component of highly 
variable nature, the singular point could not be used for 
this normalisation. The specular component in the 
transparencies not only varies from person to person and time 
to time but depends on the exposure used, since it usually is 
bright enough to saturate the film. Normalization was thus 
carried out w. r.t. an intensity derived from that measured on 
the initial, curve, which was assumed to be outside the 
specular region. 
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SHADE 298 70:84:29 
PITCH= 0.0G0 YAl/= 0.000 
DIMG= 140.0 DO8J=1000.0 


85:33:18 GOOD 2 
ROLL= 0.080 HAG= 0.7 
EYES= 0.0 



Figure 52 : Solution obtained for a nose. 

Note gaps left by the breaks caused by 
the nostrils. 



SHADE 301 


70:04:30 10:34:14 GOOD 2 


175 



Figure 54 : Four views of solution obtained for a nose. 
(With some hidden lines eliminated). 



SHADE 314 70:05:06 87:18:08 GOOD 

PITCH= 0.000 YAV/= 0.800 ROLL= 0.0G0 MAG 
DIMG= 114.0 DOBJ=l000.0' EYES= 68.0 



Figure 55 : Stereo-pair of solution obtained for a nose. 
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4 . 3.2 NORMALIZATION PROCEDURE: 

In order to simplify photographing the subjects/ it is 
necessary to make some decisions about which factors one is 
going hold fixed and which are to be taken care of by some 
normalization in the program. Altough it is possible to hold 
the head in a standard position by means of a bite-bar/ it is 
inconvenient and it Is preferable to let the program take 
care of small head-rotations. The distance from the camera 
to the subject on the other hand is very easy to determine 
and therefore no normalization of size was used. For 
pictures of the whole head such size normalization wou l d be 
fairly accurate/ whereas it cannot be for images of the nose 
alone which does not present sharp features to take 
measurements of. 

The rotational normalization procedure to be described can 
handle quite large (<TT/6) rotations in both pitch (rotation 
about an ear to ear axis) and roll (rotation about a tip-of- 
nose to back-of-head axis). Yaw (rotation about a top-of- 
head to throat axis) is restricted by the requirement that 
almost all of the surface of the nose should be visible. For 
some noses this restricts the rotation to fairly small angles 
of course this presents no problem when taking the 


photog raph . 
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LINE ALONG RIDGE 



Figure 56 : Illustration of rotational normalization 
procedure. 


i • 



Figure 57 : Illustration of parameters abstracted from one 
horizontal contour through the nose. 
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Independence of rotation is achieved by means of a routine 
which establishes the orientation of the shape calculated and 
then rotates it into a standard position. In addition the 
parameters in the final comparison procedure where chosen to 
be independent of small remaining errors in the orientations. 
The orientation of the shape calculated is estimated from two 
horizontal contours through the nose/ one passing through the 
tip of the nose, the other higher up on the ridge. These 
contours of course are only defined as sequences of points 
where the characteristics and rings pass through each plane. 

The most forward points defined by these contours are 
calculated by fitting a parabola to the three points with 
lowest z coordinates. For each of the two contours we get 
one such forward point/ connecting them we obtarn a line 
which runs approximately along the ridge of the nose. This 
line is rotated into a standard position (Lying in the y-z 
plane and leaning TT/6 from the vertical). 

The lower contour (through the tip of the nose) is also used 
to estimate rotation about the vertical axis. The two points 
on this contour at a given distance from the most forward 
point define two angles w.r.t. the z-axis. The desired 
rotation is one half of the difference of these tv/o angles. 
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The three angles so determined are small and can thus be 
treated independently. The rotation of the shape is 

performed about the center of the spherical cap used to 
determine the initial curve, i.e. a point just inside the tip 
of the nose. The whole process is repeated iteratively three 
times. The errors remaining are almost always less than 0,01 
radian (0,5 °). It was found that using only the few points 
indicated to determine the rotation was quite satisfactory, 
although better accuracy is no doubt obtainable if the 
calculation employed averages over several points. 


4,3,5 COMPARISON PROCEDURE: 

After the data has been brought into a standard orientation, 
we would like to abstract a small number of parameters which 
contain most of the information for comparison purposes. A 
rather arbitrary decision was taken to use estimates of the 
distance of the ridge of the nose from the standard line (In 
the y-z plane and leaning t t/6 from the vertical), the width 
of the nose about half-way down to the cheek and the depth of 
the cheek from the ridge of the nose. These quantities where 
measured for each of five horizontal contour planes, the 
lowest through the tip of the nose, the highest a bit below 
the saddle point (the bridge between the eyes). The fifteen 
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values so obtained are the only data used in the final 
comparison procedure. 

The distance down the side of the nose from where 
measurements of the width of the nose are taken varies with 
the contours, going from some large value for the plane 
passing through the tip of the nose to one-half that value 
for the highest contour. The distance at which the depth of 
the cheek is measured is twice that at which the width of the 
nose is measured and thus also varies from contour to 
contour. The depth of the cheek is the average of the depth 
obtained on the left side and that obtained on the right. 
The fifteen measurements obtained for each tran sparency a re 
stored in a table together with the number of the 
transparency. 

The purpose of the comparison procedure is to establish if 
any of the stored measurements match those obtained from a 
new transparency. To determine this, a pseudo-distance is 
calculated (in the 15-dimensional vector space), between each 
stored vector and the new vector. The pseudo-distance is a 
weighted r.m.s. of differences in coordinates C 7 21 , where the 
weights are proportional to the standard deviation observed 
for that coordinate. 
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d * = £ (x'-x'w 

L 

where d Is the pseudo-distance, x ? and x/ the components of 
the two vectors/ and cT. the standard deviation of the i ^ 

V# 

component. The uncertainty in the depth to the cheek is 

greater than that in the width of the nose, for example, and 
it therefore has a lower weight than the latter. This 
procedure gives a comparison test which is in some sense 
optimal [ 7 2 3 . 

No doubt other comparison procedures and other choices of 
parameters would have been equally useful; in particular it 
soon become apparent that fewer than 15 parameters would have 
been equally as selective. The point is that once one has 
data as complete as a full description of the shape, almost 
any method will work and it is not even necessary to display 
great sophistication in one's use of statistics. 


4.3.4 RESULTS OF THE NOSE-RECOGNITION PROGRAM: 

15 transparencies of 12 noses were used in this ex per iment. 
The pairs of transparencies for the three noses which were 
photographed twice differed in camera to subject distance, 
head rotation and exposure. A total of 30 shapes were 
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calculated,, two each for those noses of which only one 
transparency was available (they differed because of the 
noisy nature of the data). For each shape so determined, the 
rotational no rrra l i Eat ion was applied and the 75-tuple 
description abstracted. The pseudo-distances between all 
pairs of /5-tuples were then determined. 

The pseudo-distance between /5-tuples averaged to the 
following (the units are about 0.3 mm's r.m.s.): 

7. Between transparencies of different noses - 70. (ran^e 

2.4 - 75. 5) 

2. Between t ranspa renc i es of the same nose - 2. (range 1.4 
to 2.5) 

3. Between shapes calculated from the same transparency 

7 . (range 0.1 to 2.1) 

In all cases the distance from a given shape to a related 
shape was less than a quarter of the distance to a unrelated 
shape. Simply looking for the smallest pseudo-distance (and 
checking whether it is fairly small), thus gives an effective 
recognition procedure for this small data-set. It is clear 
that for a much larger data-set unique identification would 
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be more unlikely without improving the accuracy in the 
solution and a detailed analysis of which parameters to 
abstract for optimal recorn it ion. It would however always be 
possible to separate out some small subset of the total 
stored set of nose-descriptions with very high probability 
that the nose looked for will be in this set. Bledsoe [72] 
uses the ratio of the size of this subset to the size of the 
complete stored set as a measure of the effectiveness of the 
recognition procedure. 

Repeating the operations we described here for the other 
large frontal planes (planes with normal parallel to the z- 
axis), one would obtain a face-recogn i tion procedure. It is 
very likely that the subsets of all stored face-decr iptions 
determined by applying the above method to cheeks, chin, 
forehead and nose in turn will have only a small 
intersection. This is not to say that other information 

about the face, not obtainable from the s hape-f rom-shad i ng 
method could not add to the accuracy of such a procedure. It 
must be pointed out that some of the feature points used in 
previous attacks on the f ace-recogn i t i on problem are not 
defined by sharp discontinuities (for example the tip of the 
nose) and could best be obtained from a description of the 


shape . 
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The restriction about the positioning of the l ight-source 
could be removed if one took the trouble to measure the 
reflectivity function in more detail and either recorded the 
positioning of the light-source or worked out in detail a 
method for finding the single light-source from the shadows 
in the image (which should not be very difficult since we 
know the approximate shape of the object we are looking at). 
The full face-recognition problem was not tackled since it 
would require a great deal more work without further 
illustrating the method of determining shape-from-shading. 
Also it will be noted that the study involved a small set of 
noses - a study with a large data-set would contribute little 
more to the understanding of the method. 

Some of the difficulties encountered when determining the 
shape of noses are perhaps worth mentioning. Firstly/ most 
noses are not completely visible from any given point of 
view. Most notably the underside (between the nostrils) is 
frequently not visible, and often a small area on the side of 
the nostrils is also hidden. This forced a choice of 
parameters which did not depend on these areas, Naturally 
the information of whether these areas are visible could in 
itself be useful in the recognition procedure if it could be 
reliably determined. In fact our program does not, because 
of the combination of poor resolution in the image-dissector 
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and the simple-minded edge detector. This could be 
circumvented by placing the light-source slightly above the 
camera, thus ensuring that there always is a narrow shadow 
below the nost r i ls. 

When the solution is erroneously continued across an edge 
(such as that above the nostrils), a second undesireable 
effect appears because of the sharpening procedure. The 
incorrect coordinates of the points calculated after the edge 
is crossed have some effect on their neighbors due to this 
and thus decrease the accuracy of the solution obtained 
nea rby. 

Another problem is that some noses have not one, but two 

closely spaced tips (probably because of the underlying 

cartilage consisting of two symmetrical parts). This causes 
the characteristic growing from one of these peaks towards 
the other to stop, since it is approaching another singular 
point. A simple solution consists of chcsing the radius of 
the initial curve large enough to completely include both 
singular points. Finnally one finds that some noses 
(part icularly those belonging to females) have very low 
ridges near the eyes, making it difficult to determine a 

meaningful value for the width of the nose at that point. 
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It should be noted that the reflectivity function was not 
determined with great precision and no account was taken of 
its variation from person to person. It was not important 
that the shape calculated was very close to that of the nose 
from which the image was taken, but rather that differences 
in the shapes of noses should show up as differences in the 
calculated shapes and that shapes determined from 
transparencies of the same nose should be similar. If the 
images were all produced with the heads in the same 
rotational position, the distortions would have made no 
difference at all. For the small head rotations encountered, 
the effect of the relatively minor distortions v/as very 
sma l l . 


4.4 SUMMARY AND CONCLUSIONS: 

After defining the reflectivity function, an equation v/as 
found relating the intensity measured in the image of a 
smooth opaque object to the shape of the object. This 
equation was then shown to be a first-order non-linear 
partial differential equation in two unknowns and the 
equivalent set of five ordinary differential equations v/as 
derived. A number of especially simple cases we re 

discussed, in particular applications to lunar topography and 
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the scanning electron microscope. Methods were described for 
obtaining the auxiliary information required (e.g. the 
reflectivity function) and how to avoid the need for an 
initial known curve on the object. Of importance too is the 
method demonstrated for continuously updating p and q 
(sharpening) as the solution progresses. 

The half-dozen or so other depth-cues were ignored here to 
allow a comprehensive treatment of shading. The analytical 
approach to the problem of determining shape from shading was 
developed to demonstrate that an exact solution is possible 
and to determine just what the limitations of this approach 
are. This is not to say that a more heuristic, approximate 
approach does not have its merits too for certain types of 
objects [M3. It was decided to produce a program to allow 
experimentation with the solution method because many ideas 
in the field of artificial intelligence and visual perception 
are of little value until they can he tried on real data. 
Fortunately an image-dissector was available to provide input 
of image intensities to the computer. 

Two programs were presented, one solving the O.D.F, s for the 
characteristics sequentially, the other in parallel. 
Advantages of the latter approach were found to be several. 
Finally this latter program was adapted to provide input for 
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a n ose - recogn i t i on procedure. 

It has been made apparent that shading is valuable as a 
monocular depth-cue altough it may not be as accurate as some 
others. It must be emphasized that no claim is made that 
people employ this depth-cue in the same way. It may be that 
the human visual system does not actually determine the shape 
in three-space and if it does so it is likely that it uses a 
different method. However there will be many similarities 
between the two systems (e.g. in the errors they make) 
because they utilize the same data. 


4.4.1 SUGGESTIONS FOR FUTURE WORK: 

J. It woul d be instructive (but very tine consuming) to 
measure many reflectivity functions and see how many 
fall into the pattern of a matt component/ approximately 
varying with cos(i)/ plus a specular component. If it 
could be shown that most real reflectivity functions 
fall into this class, the method presented would be more 
useful since it could determine approximate shapes 
without knowing much more about the reflectivity 
function. 
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2. It ray be possible to find more simplifying conditions 
s.a. the ones found with certain lighting conditions, 
positions of the light-source and special reflectivity 
function s . 

3. Other solution methods may be found, or modifications to 
the integration method might increase the accuracy. 
Perhaps a difference method on a fixed grid could be 
found which somehow gets around such problems as that of 
ambiguity edges. 

4. One could study the two related problems of finding the 
reflectivity function, given the shape of the object and 
the light-source distribution and finding the light- 
source distribution given the reflectivity function and 
the shape. 

5. Further study of certain types of inconsistencies and 
their use is indicated. Here for example we find the 
problem of deciding whether certain faces in an image of 
several polyhedra could consistently belong to one 
obj ect . 

6. Some effort could be directed towards implementing more 
fully some of the ideas developed theoretically here. 
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s.a. shadow bridging, handling multiple sources and 
multiple singular points. 

7. Expanding the nose-recognition program into a full face- 
recognition program would increase its usefulness. 

%, One could study in more detail how people use the depth- 
cue of shading and how bad animals are at it. Perhaps 
one can get a better clue as to whether people develop a 
three-dimensional model of the object from the shading 
or if they use the shading information in some other 
way. 

9. There are probably a few more loose ends such as the 
problem of how to start the solution if no convex or 
concave singular points are available. Can one do 
anything at all with saddle points (even though they can 
camouflage themselves to be indistinguishable from 
simple convex or concave singular points)? 


In addition 

to interpolat ion/ 

is it 

re as onb l e 

to 

ext rapolate? 

That 

is, can 

one 

gene rate 

new 

cha racte ristics 

next to a 

so l ution 

s hee t 

to exp l ore 

new 

areas. In 

part icula r 

when a 

break 

appea rs 

In a 


solution surface can it be patched-up tater? 
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11. More methods v/ilt have to be found to deal with the 
three-dimensional structure once one has determined it. 

12. As pointed out earlier,, the use of constant size steps 

along the cha rac t e r is t ics may not be ideal (remember 
that we can adjust the step-size by choosing a different 
-X ). One particularly attractive idea would be to use 
steps corresponding to constant Intensity change in the 
image. This would turn the rings into constant 

intensity contours, rather than curves of constant arc- 
distance from the singular point. 

13. Many objects have surfaces whose reflectivity cannot be 
described by a function of three angles, or are so 
specular that our methods are of little avail. One 
might try to discover methods of dealing with such 
objects. Examples are chrome ca»—bumpers, translucent 
wax, hair and a glass of water. 



Po^e 195 


5. REFERENCES: 

J. V. P. Fesenkcv: 'PHOTOMETRIC INVESTIGATIONS OF THE 

LUNAR SURFACE' 1929 Astronomochheskii Zhurnal 5 
(Translated by Redstone Scientific Information Centre). 

2. J. van Diggelen: 'A PHOTOMETRIC INVESTIGATION OF THE 
SLOPES AND HEIGHTS OF THE RANGES OF HILLS IN THE MARIA 
OF THE MOON' July 1951 - Bulletin of the Astronomical 
Institute of the Netherlands, 

3. D. E. Willingham: 'THE LUNAR REFLECTIVITY MODEL FOR 

RANGER BLOCK III ANALYSIS' November 1964 - Technical 

Report 3 2-664 Jet Propulsion Laboratory. 

4. T. Rindfleisch: 'PHOTOMETRIC METHOD FOR LUNAR 

TOPOGRAPHY' March 1 966 - Phot op; rammet r ic En£ inee r Inn , 

5. D. R. Garabedian: 'PARTIAL DIFFERENTIAL EQUATIONS' 1964 
- John Wiley. 

6. R. D. Richt me ye r and K. W. Morton: 'DIFFERENCE METHODS 
FOR INITIAL VALUE PROBLEMS' 1 957 - John Wiley 


", » » » •» ■ v ■ ' >' I 


. I Pagp 1 9 6 

7. R. W. Hamming: 'NUMERICAL METHODS FOR SCIFMTISTS AMD 
ENGINEERS' 7962 - McGrav; Hill. 

8. K. C. Kelly: 'A COMPUTER GRAPHICS PROGRAM FOR THE 

GENERATION OF HALF-TONE IMAGES WITH SHADOWS'. November 
1 9 69 - University of Illinois. 

9. B. K. P. Horn: 'THE IMAGE DISSECTOR 'EYES" August 1969 
- A.I. Memo 17 8 - M.I.T. 

10. S. A. Coons: 'SURFACES FOR COMPUTER-AIDED DESIGN OF 

SPACE FIGURES' February 1968 - Mechanical Eng M.I.T. 

11. P. R. Thornton: 'THE SCANNING ELECTRON MICROSCOPE' 

November 1965 - Science Journal 

12. W. W. Bledsoe: 'MAN-MACHINE FACIAL RECOGNITION' August 
1 966 - Panoramic Research Report 22 - Palo Alto 

13. C. Engelman: 'MATHLAB' August V968 - IFIP Congress 

[page B97 to B963 - Edinborough 

14. L. J. Krakauer: 'COMPUTER ANALYSIS OF VISUAL PROPERTIES 
OF CURVED OBJECTS' June 1970 - Electrical Eng M.I.T. 



